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INVESTIGATING  SORBENT  GEOMETRY  RESPONSE 
SURFACES  FROM  SORPTION  RATE  DATA 

ABSTRACT 

AFIT/GO  A/EN  V/98M-0 1 

Sorption  and  desorption  of  hydrophobic  organic  compounds  (HOCs)  into  soil 
particles  occurs  at  Air  Force  contamination  sites.  Long-term  desorption  extends  clean  up 
time,  costing  billions.  Accurately  modeling  desorption  will  reduce  costs  and  improve  clean 
up  designs.  State  of  the  art  models  depict  soil  as  uniform  spherical  particles  which  lose  the 
effect  of  longer  sorption  path  lengths.  An  alternate  approach,  the  multiple  sites  in  series 
(MSS)  model,  describes  the  sorption  capacity  of  a  mixture  of  soil  particle  sizes  and  shapes 
using  a  composite  particle  defined  by  a  two  parameter  statistical  distribution.  When  the  MSS 
model  uses  a  general  radial  geometry  function  as  the  statistical  distribution,  unique 
parameters  can  not  be  fit  to  laboratory  experimental  data.  This  research  employs  response 
surface  methodology  (RSM)  to  investigate  the  nonunique  parameter  model  deficiency. 
Analysis  shows  that  the  response  surface  for  the  general  radial  geometry  function  is  a  trough 
with  no  unique  minimum.  A  unique  minimum  is  found  when  the  statistical  distribution  is 
changed  to  a  gamma  distribution. 
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I.  INTRODUCTION 


Hydrophobic  organic  compounds  contaminate  groundwater  at  many  hazardous  waste 
sites.  These  chemicals  threaten  the  health  and  safety  of  those  where  the  groundwater  is  used 
for  water  supply.  A  great  deal  of  time,  money,  and  research  is  being  applied  to  clean  the 
chemicals  from  the  groundwater,  but  the  efforts  are  not  always  successful. 

To  make  remediation  efforts  more  successful,  research  is  focusing  on  the  modeling  of 
solute  transport.  The  transport  process  can  be  classified  into  three  broad  categories:  1) 
solute  transfer  between  phases  (sorption,  dissolution),  2)  spatial  transport  of  a  solute  within  a 
fluid  (advection,  dispersion,  diffusion),  and  3)  chemical  change  of  solutes  (biodegradation). 
This  research  effort  examines  solute  transfer  between  phases. 

The  purpose  of  a  model  simulating  solute  transport  is  to  predict  the  contaminant 
concentration  at  a  particular  time  and  point  in  space  (Konikow,  1981).  However,  there  are  so 
many  parameters  to  estimate  that  the  task  is  difficult  for  one  model  to  be  able  to  compute  the 
values  effectively.  For  this  and  other  reasons,  several  models  are  used  to  calculate  parameter 
estimates  that  are  used  in  the  overall  calculation  of  a  contaminant  concentration. 

A  chemical  contaminant  in  water  will  migrate  toward  locations  where  it  is  less 
concentrated.  This  process  is  known  as  diffusion.  In  soil  systems  where  the  concentrations 
change  with  time,  Fick’s  second  law  is  applied.  In  one  dimension,  Fick’s  second  law  is 
written  as 

dC  d2  C 

d  t  d  x2 

where  C  =  concentration, 
t  =  time,  and 
D  =  diffusion  coefficient. 
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In  soil  systems,  diffusion  will  not  proceed  as  fast  as  in  water  because  the  solute  must  travel 
longer  distances  due  to  the  tortuous  paths  around  solids.  To  account  for  this,  an  effective 
diffusion  coefficient,  Deff,  is  used,  which  includes  a  measure  of  the  effect  of  the  shape  of  the 
flowpath  followed  by  the  molecules  (Fetter,  1992).  Calculation  of  the  effective  diffusion 
coefficient  is  simplified  by  using  data  obtained  in  laboratory  experiments.  Computer  models 
can  fit  Deff  coefficients  using  known  concentrations  from  experimental  data  over  time. 

Heyse’s  (1994)  multiple  sites  in  series  (MSS)  model  will  fit  several  parameters  to 
experimental  data,  one  of  which  is  Deff.  The  MSS  model  can  also  fit  geometry  parameters  to 
data  sets,  which  allows  characterization  (size,  shape)  of  a  composite  soil  particle.  Previous 
research  has  found  two  problems  with  the  MSS  model,  namely,  that  the  Deff  parameter 
depends  on  the  size  of  the  synthetic  sorbent  spheres  used  in  laboratory  experiments,  and  the 
geometry  parameters  found  for  mixtures  of  particle  sizes  are  not  unique. 

In  order  for  modeling  to  be  useful,  consistent  results  should  be  obtained.  The  size 
dependent  Deff  coefficients  are  unexpected.  For  a  given  system,  the  diffusion  coefficient 
should  generally  not  be  dependent  on  the  size  of  the  particles  for  correctly  defined 
geometries.  An  error  is  occurring  that  causes  the  size  dependent  Deff  coefficients.  This  error 
could  involve  the  assumption  of  Fickian  diffusion,  the  laboratory  data,  the  composite  shape 
approach,  or  even  the  model  code  itself,  just  to  name  a  few  of  the  possible  sources.  If  the 
MSS  model  is  to  be  used  to  help  predict  contaminant  concentrations,  it  must  also  be  able  to 
predict  unique  geometry  parameters  for  a  given  soil  system  so  that  different  systems  can  be 
characterized  in  laboratory  experiments. 

To  improve  the  MSS  model,  this  research  will  examine  and  explain  the  size 
dependence  of  the  Deff  coefficients.  The  utility  of  the  model  will  increase  if  the  diffusion 
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coefficients  it  predicts  are  within  a  range  of  expected  values  and  also  constant  for  a  given 


sorbent. 

The  second  objective  of  this  research  is  to  characterize  the  nature  of  the  geometry 
parameters’  response  surface.  The  sum  of  squares  for  error  (SSE)  is  calculated  between 
known  data  and  model  predictions  by  the  MSS  model.  Response  surface  methodology  will 
be  used  to  show  the  shape  of  the  SSE  response  surface  formed  by  the  two  input  geometry 
parameters.  The  shape  of  the  surface  may  explain  why  there  is  no  unique  best  fit  pair  of 
parameters.  Finally,  this  effort  will  investigate  another  function  in  the  MSS  model  to 
construct  the  SSE  surface  and  determine  if  a  unique  pair  of  geometry  parameters  can  be 
found. 
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II.  LITERATURE  REVIEW 


The  most  common  hazardous  contaminants  at  Air  Force  waste  sites  are  hydrophobic 
organic  compounds  (HOCs),  such  as  chlorinated  solvents  and  aircraft  fuel  components. 

Current  research  is  concentrating  on  how  to  model  the  movement  of  these  chemicals  in  order 
to  determine  the  best  remediation  measures.  For  this  research  to  be  effective,  knowledge  of 
the  physical,  chemical,  and  biological  processes  occurring  to  the  contaminants  in 
groundwater  is  essential. 

Sorption  Processes 

Sorption  is  an  important  physical-chemical  process  and  is  the  means  by  which 
chemicals  associate  with  solid  phases  (Schwarzenbach  et  al. ,  1993).  The  generic  term 
sorption  includes  both  adsorption  and  absorption.  Adsorption  is  the  accumulation  occurring 
at  an  interface  while  absorption  is  the  partitioning  between  two  phases.  For  absorption,  the 
equilibrium  partition  coefficient  (Kp)  quantifies  the  solute  distribution  between  the  liquid  and 
solid  phases.  If  a  soil  has  a  significant  natural  organic  carbon  content  (>  0.1  percent), 
sorption  by  the  soil  organic  matter  (SOM)  dominates  all  other  HOC  sorption  processes 
(Karickhoff,  1981). 

Sorption  is  a  complex  process  that  is  receiving  much  attention  in  the  literature.  Better 
understanding  of  this  process  will  allow  improved  model  construction  and  remediation 
design,  as  well  as  improved  knowledge  of  waste  site  contaminant  transport.  Currently, 
sorption  is  modeled  with  equilibrium  or  nonequilibrium  models,  each  with  their  own 
assumptions  and  limitations. 

A  common  assumption  made  in  modeling  sorption  is  that  of  local  equilibrium.  The 
local  equilibrium  assumption  (LEA)  states  that  the  solute  is  at  equilibrium  between  all  phases 
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at  any  point  in  space  and  is  used  to  simplify  the  transport  equations  (Valocchi,  1986;  Goltz 
and  Oxley,  1991).  According  to  the  LEA,  the  process  of  a  contaminant  partitioning  between 
sorbed  and  dissolved  phases  is  instantaneous.  However,  there  are  numerous  cases  where 
equilibrium  conditions  do  not  prevail  in  the  subsurface.  Nonequilibrium  sorption  of  HOCs  at 
the  grain  scale  exists  because  mass  transfer  between  moving  groundwater  and  SOM  is 
governed  by  diffusion,  a  slow  process  compared  to  advective  transport  with  groundwater 
flow  (Brusseau  and  Rao,  1989;  Brusseau  et  al.,  1991). 

The  two  most  accepted  theories  of  HOC  sorption  nonequilibrium  are  intraorganic 
matter  diffusion  (IOMD)  and  retarded  intraparticle  diffusion  (RIPD).  IOMD  involves 
absorption  (partitioning)  of  a  chemical  species  within  the  matrix  of  SOM  (Chiou  et  al.,  1979; 
Chiou  etal.,  1983).  RIPD  involves  aqueous-phase  diffusion  and  sorption  of  contaminant 
species  within  the  microscopic  pores  of  mineral  grains  (Wu  and  Gschwend,  1986;  Ball  and 
Roberts,  1991b). 

Nonequilibrium  sorption  is  important  because  it  results  in  both  earlier  arrival  and 
longer  tailing  of  solutes  moving  past  a  location  in  space.  There  are  several  important 
implications  as  sorption  mass  transfer  slows:  1)  the  duration  of  exposure  to  a  solute 
increases,  potentially  increasing  health  risks,  2)  more  time  is  required  for  engineered 
restoration  activities,  increasing  clean  up  costs,  and  3)  solutes  spend  more  time  in  the  sorbed 
phase  where  they  are  probably  not  available  for  biodegradation  (Scow  and  Alexander,  1992), 
possibly  decreasing  the  effectiveness  of  natural  attenuation  or  intrinsic  remediation. 
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Sorption  Models 

Sorption  kinetics  are  modeled  in  one  of  three  ways,  depending  on  the  assumptions 
made  and  the  level  of  mathematical  detail  desired.  If  equilibrium  can  be  assumed,  the  rate  of 
change  of  sorbed  phase  contaminant  is  proportional  to  the  rate  of  change  of  aqueous  phase 
contaminant,  shown  in  equation  (1) 


dS 

dt 


(1) 


where  S  =  sorbed  phase  concentration, 

Kp  =  equilibrium  partition  coefficient,  and 
C  =  solute  concentration. 

Nonequilibrium  models  employ  either  first-  or  second-order  approximations  of 
diffusion  (Parker  and  Valocchi,  1986;  Wu  and  Gschwend,  1986;  Brusseau  and  Rao,  1989; 
Ball  and  Roberts,  1991b).  The  first-order  model  uses  a  first-order  mass  transfer  expression 
to  describe  the  sorption  domain  as  a  discrete  distribution.  This  model  is  described 
mathematically  using  equation  (2) 

f=k;(KpC-S)  (2) 


where  =  first-order  mass  transfer  rate. 

The  first-order  description  of  the  sorption  domain  is  highly  simplified  which  allows 
computational  ease,  but  makes  correlations  for  new  combinations  of  chemicals  and  solids 
from  previously  fit  combinations  difficult  (Wu  and  Gschwend,  1986).  The  second-order 
description  is  based  on  molecular  (Fickian)  diffusion  and  phase  partitioning.  The  time  rate 
of  change  of  sorbed  compound  per  unit  volume  is  expressed  using  Fick’s  second  law  of 
diffusion  in  radial  coordinates,  most  often  employing  a  spherical  geometry.  This  radial 
diffusion  is  described  by  equations  (3)  through  (7) 
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(3) 


as  es  acs 
at _  p  at 


Cs=— r-fRr2Csdr 
s  ^3  Jo  s 


acs  _  Deff  a 
at  r2  dr 


3C, 

ar 


Cs(r 

ac^ 

ar 


=  R)  =  Kp 


=  0 


r=0 


c 


(4) 

(5) 

(6) 

(7) 


where  Deg  =  effective  diffusion  coefficient, 

C5  =  volumetric  concentration  in  soil  organic  matter  at  position  r  and  time  t, 

Cs  =  average  sorbed  phase  concentration, 

R  =  radius  of  sphere, 
p  =  soil  bulk  density, 

9S  =  volumetric  water  content. 

As  with  all  mathematical  modeling,  these  equations  contain  limitations.  The  first- 
order  model  requires  the  assumptions  that  1)  all  sorption  sites  are  located  at  the  same  path 
length,  and  2)  concentration  gradients  are  linear  between  the  sorbed  and  aqueous  phases. 
This  is  not  always  the  case.  The  first-order  model  does  a  poor  job  at  reproducing  long-term 
slow  desorption,  and  mass  transfer  rate  coefficients  have  been  experimentally  found  to 
depend  on  groundwater  velocity.  Also,  first-order  models  usually  assume  a  single  mass 
transfer  rate  coefficient  for  the  sorption  and  desorption  of  contaminants  to  and  from  natural 
sediments  (Brusseau  and  Rao,  1989).  This  assumption  often  causes  model  failure  when 
simulating  sorptive  mass  transfer  under  changing  solvent  velocities  (Brusseau,  1992;  Heyse, 
1994;  Heyse  et  al.,  1997)  or  extended  contamination  contact  times  (Karickhoff  and  Morris, 
1985;  Connaughton  et  al.,  1993;  Young  and  Ball,  1995;  Culver  etal.,  1997).  First-order 
models  can  not  deal  with  heterogeneity  either.  Some  research  has  accounted  for 
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heterogeneity  by  using  multiple  first-order  sites  (Connaughton  et  al.,  1993;  Heyse,  1994; 
Culver  et  al.,  1997;  Haggerty  and  Gorelick,  1995).  Although  the  multiple-particle  class 
formulation  of  the  first-order  mass  transfer  model  allows  for  variations  in  sorption  rate  and 
equilibrium,  measuring  these  parameters  is  not  experimentally  possible  (Pedit  and  Miller, 
1994).  The  first-order  models  are  approximating  diffusion,  but  a  second-order  model  will 
always  approximate  diffusion  more  closely. 

The  second-order  model  presented  above  assumes  that  all  sorption  sites  are 
distributed  as  uniform  spheres.  This  is  certainly  not  the  case  because  natural  sorbents  are 
heterogeneous  mixtures  with  complex  geometries  (Brewer,  1964;  Weber  et  al.,  1992).  Also, 
since  the  diffusion  model  requires  the  calculation  of  solute  concentrations  at  all  radial 
coordinates  in  the  sphere,  it  is  mathematically  more  complex  than  the  first-order  approach. 
Dealing  with  Heterogeneity 

The  assumption  of  homogeneous  spherical  particle  geometries  as  usually  employed  in 
second-order  diffusion  based  models  is  unrealistic.  SOM  is  very  heterogeneous,  being 
comprised  of  various  molecular  sizes  and  densities.  Schnitzer  (1978)  and  Stevenson  (1982) 
report  SOM  to  be  a  flexible,  cross-linked,  branched,  amorphous,  polyelectrolytic  polymeric 
substance  with  various  distributions  of  grain  sizes.  A  soil’s  diagenetic  history  involves  the 
physical  and  chemical  changes  occurring  during  and  after  formation.  Weber  et  al.  (1992) 
and  Young  and  Weber  (1995)  report  significant  differences  in  soil  based  on  diagenetic  age  of 
the  soil,  and  hence,  the  SOM  contained  in  it.  Older  SOM  has  a  condensed,  almost  crystalline 
structure  that  is  significantly  more  hydrophobic,  causing  slower  diffusion.  On  the  other 
hand,  younger  SOM  has  a  rubbery,  more  open  structure  that  is  less  hydrophobic  and  allows 
for  faster  sorption.  Clearly  SOM  is  not  homogeneous. 
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To  improve  sorption  mass  transfer  models,  the  effects  of  heterogeneity  are  usually 
included  by  using  one  of  three  methods.  The  first  technique  is  to  model  all  soil  particles  as 
an  equivalent  sphere.  The  equivalent  sphere  can  be  determined  in  one  of  two  ways.  The  first 
is  to  find  the  equivalent  sphere  that  gives  the  same  specific  interfacial  area  as  the  aggregated 
sorbent.  The  second  is  to  select  a  sphere  of  mean  radius  for  the  soil  under  consideration  and 
fit  a  diffusion  coefficient.  The  mean  radius  is  often  selected  as  the  sieve  size  through  which 
fifty  percent  of  the  mass  will  pass.  The  spherical  diffusion  model  can  then  be  fit  to  data  by 
modifying  the  diffusion  coefficient.  However,  by  describing  concentrations  at  all  path 
lengths  in  the  sphere,  all  models  require  a  significant  amount  of  computational  effort  to  solve 
(requires  larger,  slower  models).  Another  problem  that  arises  is  that  the  effect  of  longer  path 
lengths  are  lost  (Haggerty  and  Gorelick,  1995).  Likewise,  the  multiple-particle  class 
formulation  of  the  Fickian  diffusion  model  can  only  address  variability  down  to  the 
individual  particle  scale,  which  assumes  an  individual  particle  to  be  homogeneous  and  does 
not  allow  for  intraparticle  variations  in  sorption  rate  and  equilibrium  (Pedit  and  Miller, 

1994).  These  models  are  also  very  complex. 

The  second  method  used  to  account  for  heterogeneity  involves  using  the  continuum 
of  compartments  concept  proposed  by  Connaughton  et  al.  (1993).  The  multiple  sites  in 
parallel  (MSP)  model,  as  shown  in  Figure  1,  employs  multiple  compartments.  Each 
compartment  is  defined  by  a  first-order  mass  transfer  rate.  Connaughton  et  al.  (1993)  and 
Culver  et  al.  (1997)  use  a  gamma  distribution  to  describe  the  mass  transfer  rates,  while 
Heyse  (1994)  and  Culver  et  al.  (1997)  use  a  lognormal  distribution,  all  with  varying  degrees 
of  success.  MSP  models  of  this  type  are  not  computationally  difficult,  but  analytical 
solutions  are  only  possible  for  the  simplest  of  boundary  conditions.  The  biggest  problem 
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found  is  that  changing  the  boundary  conditions  results  in  different  distributions  (model 
parameters)  (Connaughton  etal.,  1993;  Heyse,  1994). 


0 

Figure  1,  Multiple  Sites  in  Parallel  (MSP)  Model  (Heyse,  1994) 

A  variation  of  the  MSP  model  was  described  by  Haggerty  and  Gorelick  (1995).  Their 
multirate  model  is  different  from  others  in  that  instead  of  a  statistical  distribution  of  rate 
coefficients,  the  size  and  rate  constants  are  derived  so  that  they  give  results  identical  to 
diffusion-based  models.  Multiple  sites  or  compartments  each  have  their  own  first-order  rate. 
The  rates  in  each  compartment  can  be  fixed  so  that  the  model  gives  the  same  results  as  a 
layered,  cylindrical,  or  spherical  diffusion  model.  Different  particle  sizes  and  shapes  can  be 
modeled  simultaneously  by  weighting  each  particle  size/shape  by  its  fraction  of  sorbing 
capacity.  A  potential  drawback  is  that  the  solution  is  derived  using  a  uniform  concentration 
as  the  initial  condition  in  the  sorbed  phase.  This  initial  condition  may  result  in  different 
solutions  when  the  boundary  conditions  change. 
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The  third  method  of  dealing  with  heterogeneity  in  diffusion  models  is  to  include  all 
geometric  shape  descriptions  in  one  composite  (but  nonspherical)  shape.  Heyse  (1994) 
introduced  the  multiple  sites  in  series  (MSS)  model  as  depicted  in  Figure  2.  It  is  simply  a 
finite  difference  solution  to  the  diffusion  equation  into  any  geometry  as  described 


0  8 


Figure  2,  Multiple  Sites  in  Series  (MSS)  Model  (Heyse,  1994) 


by  a  frequency  distribution  of  sorption  capacities  along  a  single  diffusion  path,  de  Venoge 
(1996)  and  Mika  (1997)  used  the  geometric  function  given  in  equation  (8)  to  describe  a 
composite  particle’s  geometry.  By  including  the  shape  factor  (A.)  in  the  frequency 
distribution,  different  particle  geometries  can  be  described. 


f  (8)  = 


A.+1 


>.+  t 


(R-8)‘ 


0<8<  R 
otherwise 


(8) 


where  f(8)  =  frequency  distribution  representing  sorption  capacity, 

A  =  shape  factor, 

R  =  maximum  diffusion  path  length,  and 
8  =  variable  diffusion  path  length. 

A  uniform  coating  has  X  =  0,  X  =  1  describes  a  cylinder,  and  a  spherical  particle  results  from 
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X  =  2.  Equation  (8)  represents  the  cross-sectional  area  normalized  to  the  volume  of  the 
sorbent  over  the  range  from  6  =  0  (sorbent  /  water  interface)  to  8  =  R  (center  of  the  particle). 
Note  that  the  diffusion  path  length  8  is  measured  from  the  liquid  interface  in  equation  (8),  but 
equations  (3)  -  (7)  measure  radius  from  the  center  of  the  particle.  Equation  (8)  makes  the 
diffusion  path  length  computationally  easier  to  use.  It  describes  a  distribution  of  sorption 
sites  along  a  diffusion  path  length.  Additionally,  higher  order  geometries  can  be  represented 
by  increasing  the  shape  factor.  Higher  shape  factors  result  in  a  majority  of  sorption  sites  at 
short  diffusion  path  lengths,  but  with  some  sites  at  very  long  path  lengths.  The  description  of 
diffusion  path  lengths  with  a  frequency  distribution  that  includes  a  shape  factor  is  intuitively 
appealing  since  it  can  account  for  the  complex  geometries  associated  with  natural  solid 
phases.  However,  as  shown  by  Mika  (1997),  the  disadvantage  of  this  distribution  is  that  for 
shape  factors  (A.)  much  greater  than  two,  the  distribution  gives  a  nonunique  solution. 
Desirable  Distribution  Characteristics 

To  expand  on  the  work  of  de  Venoge  (1996)  and  Mika  (1997),  this  thesis  tests  the 
gamma  distribution  for  use  in  the  MSS  model.  For  a  statistical  distribution  to  be  desirable, 
there  should  be  no  more  than  two  fitting  parameters.  This  prevents  the  model  from 
becoming  too  complex.  These  fitting  parameters  need  to  describe  the  composite  geometry  of 
soil  particles  using  an  f  (8)  distribution.  The  shape  of  the  distribution  must  provide  high 
frequencies  at  low  8  values  to  account  for  a  rapid  uptake  of  solute  and  a  long  tail  to  describe 
rate-limited  sorption  over  extended  time  periods. 

The  distribution  must  also  be  unique  for  a  given  sorbent.  The  distribution  evaluated 
at  zero  is  the  specific  interfacial  area,  as  defined  below 
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f  (0)  =  a  = 


interfacial  area 


volume 

The  first  normalized  moment  (mean)  of  the  distribution  is  the  average  diffusion  path  length, 


[°°8-f(8)d8 

8  where  8  =  — - =  f°°  8  •  f  (8)  d8 ,  since  the  denominator  is  one.  Mika  (1997) 

’  1*00  JO 

J0  f(8)dS 


found  that  for  any  best  fit  geometry, 


the  value  D=  was  constant. 
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Candidate  Distributions 

Several  statistical  distributions  meet  the  above  criteria.  This  research  effort  tests  the 
gamma  distribution  for  use  in  the  MSS  model  and  evaluates  this  distribution  for  uniqueness 
of  fit  using  response  surface  methodology  (RSM).  See  Hastings  and  Peacock  (1975),  or  Law 
and  Kelton  (1991)  for  more  detailed  discussions  of  statistical  distributions.  The  lognormal 
and  Weibull  density  functions  are  included  here  because  they  represent  the  general  shape  of  a 
desired  distribution.  Although  this  thesis  does  not  test  their  use  in  the  MSS  model,  they 
could  have  computational  advantages  and  should  be  the  subject  of  future  research. 

The  gamma  distribution  is  depicted  below  in  Figure  3,  where  gamma  (x,  a)  is  the 
form  of  the  curves  that  are  shown.  Figure  3  only  displays  changes  in  the  shape  parameter,  a, 
which  demonstrates  the  variability  of  the  gamma  distribution.  It  is  characterized  by  a  shape 
parameter,  a,  and  a  scale  parameter,  p,  both  of  which  must  be  greater  than  zero.  As  is  seen 
in  the  figure,  when  the  shape  parameter  a  is  less  than  one,  the  density  takes  on  infinite  values 
as  x  approaches  zero.  For  certain  types  of  soil  geometries,  this  may  prove  beneficial,  since 
the  mean  of  the  diffusion  path  lengths  may  be  located  at  various  values  of  8.  As  is  shown  in 
the  figure,  the  gamma  distribution  can  take  on  many  shapes,  depending  on  the  values  of  a 
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and  p;  The  exponential  distribution  is  a  special  case  of  the  gamma  distribution,  as  is  seen  by 
gamma  (x,  1)  in  Figure  3 

gamma  (x,  0.5) 
gamma  (x,  1) 
gamma  (x,  1.5) 
gamma  (x,  3) 

x 

Figure  3,  Gamma  Probability  Density  Function 


where  F  (a)  =  J~e_u  ua_1  du ,  the  mean  is  given  by  ap  and  the  variance  by  ct-P2. 

Figure  4  shows  the  lognormal  distribution.  The  curves  on  the  graph  are  of  the  form 
lognorm  (x,  |i,  a).  This  distribution  is  also  characterized  by  a  shape  parameter,  o,  and  a  scale 
parameter,  jx,  where  o  >  0  and  -°°  <  fi  <  oo.  Notice  in  Figure  4  that  lognormal  curves  can  be 
shifted  away  from  the  origin.  As  with  the  gamma  distribution,  this  may  prove  beneficial  for 
certain  types  of  soil  geometries. 
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The  lognormal  density  function  is  given  in  equation  (10) 


f(x) 


^2ko‘ 


exp 


-(lnx-|i)' 

2o2 


(10) 


with  mean  (i-exp(0.5-a2)  and  variance  |i2-(o2-l)-exp(<r). 


Figure  4,  Lognormal  Probability  Density  Function 


The  last  distribution  to  be  considered  is  the  Weibull,  shown  in  Figure  5.  Graphed 
curves  are  of  the  form  Weibull  (x,  a)  where  only  changes  in  the  shape  parameter,  a,  are 
demonstrated. 

As  with  the  gamma  distribution,  if  the  shape  parameter  is  less  than  one,  the  density 
takes  on  infinite  values  as  x  approaches  zero.  Certain  types  of  soil  geometries  may  require 
this  capability.  Also  note  that  as  the  shape  parameter  increases,  the  Weibull  distribution 
begins  to  resemble  a  normal  distribution.  The  Weibull  distribution  has  shape  and  scale 
parameters,  a  and  P  respectively,  both  greater  than  zero.  The  Weibull  density  function  is 
given  by  equation  (11) 
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Figure  5,  Weibull  Probability  Density  Function 
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III.  METHODOLOGY 


Data  Sets 

This  research  compares  the  ability  of  Heyse’s  (1994)  MSS  model  to  describe  sorption 
by  heterogeneous  mixes  of  sorbents  as  representative,  composite  particles  using  two  different 
functions.  The  general  radial  geometry  function  (equation  8)  and  the  gamma  density 
function  (equation  9)  are  used  to  describe  the  sorbent  capacity  of  a  representative  composite 
particle  along  a  diffusion  path  length.  These  models  are  tested  using  two  data  sets. 

One  data  set  is  a  simulation  of  sorptive  uptake  of  tetrachloroethene  (PCE)  by  a  sandy 
aquifer  material  from  Borden,  ON  in  a  batch  system.  This  simulation  is  created  using  the 
model  of  Haggerty  and  Gorelick  (1995)  and  based  on  Borden  soil  characteristics  as  defined 
by  Ball  and  Roberts  (1991a,  b),  shown  in  Table  1.  The  data  set  produced  from  this 
simulation  is  presented  in  Appendix  A.  It  is  used  as  input  for  the  MSS  model  as  a  natural 
sorbent. 


Table  1,  Borden  Site  Data  (adapted  from  Ball  and  Roberts,  1991a,  b) 


Particle  Fraction 

Mass  Fraction,  F 

R,  mm 

Kd,  mL/g 

Ds,  crrtVsecxlO1 1 

1 

0.0091 

0.6 

8 

11.2 

2 

0.0524 

0.3 

2.9 

8.28 

3 

0.163 

0.16 

1.2 

5.89 

4 

0.257 

0.11 

0.55 

3.27 

5 

0.315 

0.075 

0.3 

5.29 

6 

0.615 

0.048 

0.26 

3.92 

7 

0.0341 

0.03 

0.82 

1.26 

Composite 

0.9956 

0.727 

3.39 

Bulk 

1 

0.11 

0.74 

3.39 
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The  second  data  set  consists  of  sorption  of  anthracene  in  a  methanol/water  solution 
by  paraffin  beads  conducted  by  Mika  (1997).  The  manufacture  of  the  paraffin  beads,  a 
synthetic  sorbent  created  for  the  study  of  multisite  sorption  models,  is  outlined  by  de  Venoge 
(1996).  The  size  and  consistency  of  these  paraffin  beads  is  shown  in  Table  2.  Mika  (1997) 
uses  these  same  beads  to  conduct  sorption/desorption  experiments  in  a  batch  system  on  two 
different  mixtures  of  paraffin  bead  sizes.  The  even  mix  contained  7  large,  12  medium,  and 
23  small  beads.  The  uneven  mix  contained  2  large,  5  medium,  and  53  small  beads.  The 
equilibrium  partition  coefficient  (Kp)  is  estimated  to  be  9.35  mL/g  by  linear  regression 
(R2  =  0.9997).  Mika’s  (1997)  even  data  set  is  included  in  Appendix  B. 


Table  2,  Synthetic  Soil  Statistics  (de  Venoge,  1996) 


Statistics 

Sphere  Size 

Large 

Medium 

Small 

Average  Mass  (g) 

0.09253 

0.05511 

0.02808 

Standard  Deviation 

0.00342 

0.00158 

0.00110 

Maximum 

0.09740 

0.05972 

0.03010 

Minimum 

0.08610 

0.05232 

0.02491 

Count 

21 

36 

70 

Max  Radius  (cm) 

0.297 

0.249 

0.1997 

Deff  (cm2/s  x  10'8) 

2.12 

2.01 

1.96 

Composite  Sorbent  Particle 

Natural  sorbent  particles  are  heterogeneous;  the  approach  of  this  research  is  to 
incorporate  heterogeneity  into  a  composite  sorbent  particle.  This  is  illustrated  for  two 
separate  cases:  1)  homogeneous  material  (same  Kp,  Deff),  heterogeneous  shape,  and  2) 
heterogeneous  material  and  shape. 

Case  1:  Homogenous  material,  heterogeneous  shape.  Assume  each  particle  class 
(size  and  shape)  can  be  defined  by  its  mass  fraction,  Fj,  its  radius,  Rj,  and  its  shape  factor,  Xj. 
Since  each  mass  fraction  has  the  same  Kp  and  Deff,  the  two  distributions  can  be  added  to  form 
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one  composite  sorbent  particle  distribution,  as  shown  in  Figure  6,  where  fj  (A,, 8)  is  as  defined 
in  equation  (8). 


- Fjfj  (0, 8) 

- Fjfj  (1,8) 

- Fj-fj  (2, 8) 

- fc(A,8) 


Figure  6,  Particle  Distributions 


The  distribution  of  the  composite  particle  fc(8)  is  defined  in  equation  (12) 

fc(8)=XjFJ-fj(8)  (12) 

Once  the  composite  distribution  is  determined,  a  composite  particle  can  be  obtained  by  using 
the  equations  for  specific  interfacial  area  (a)  and  average  diffusion  path  length  (s)  from 
Chapter  II. 

Case  2:  Heterogeneous  material  and  shape.  The  sorption  characteristics  (Kp,  Deff)  of 
natural  sorbents  will  vary  along  with  particle  size  and  shape.  Individual  particles  must  be 
normalized  to  an  average  Kp  and  representative  Deff  so  a  single  composite  particle  can  be 
defined.  The  average  Kp  is  obtained  by 
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(13a) 


K 


p,  avg 


(13b) 


Several  variables  need  to  be  transformed.  A  prime  symbol  will  denote  a  virtual 
particle.  The  virtual  mass  and  mass  fractions  are  defined  in  equations  (15)  and  (16). 


K, 


Mj  =>Mj' 


P.J 


Ff 


Kp.avg 

Kp,j 


K, 


p.avg 


In  order  for  the  mass  flux  of  the  real  and  virtual  particles  to  remain  the  same 

Deff.j 


a 

a 

X 

II 

O 

Q 

X 

(14) 


(15) 


(16) 


where  Do  is  a  representative  diffusion  coefficient. 

The  specific  interfacial  area  (a)  and  the  average  diffusion  path  length  (b)  are  defined 


by 


X  i  + 1 

aj  =  fj(0)=-^- 


R: 


(5)d6=ri 


X  j+ 2 


(17) 


(18) 


The  particle  shape  must  be  maintained  in  the  virtual  particle,  so  X  i  =  X  j .  By  inserting 
equations  (17)  and  (18)  into  equation  (16)  and  simplifying,  the  virtual  radius,  Rj,  is  shown 
in  equation  (19) 
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RrRrJ%^ 


(19) 


The  virtual  particles  are  illustrated  in  Figure  7.  Since  Kp,  avg  and  D0  are  the  same  for  each 
particle,  a  composite  sorbent  particle  distribution,  fc(8),  is  obtained  using  equation  (12)  as  in 
case  1 . 


Real 

Virtual 

Real 

Virtual 

o 

*  \ 

*  \ 

*  * 

t  • 

%  i 

*  * 

O 

n 

Particle  Class  1 

Particle  Class  2 

Kp, 

Kp,  avg 

Kp2 

Kpavg 

Deff 

Do 

Deff 

Do 

R, 

Ri 

R2 

R'2 

F, 

F, 

f2 

f2 

Figure  7,  Composite  Sorbent  Particle,  Different  Kp  and  Deff 


Response  Surface  Methodology 

RSM  comprises  mathematical  and  statistical  techniques  to  develop,  improve,  or 
optimize  various  processes  (Myers  and  Montgomery,  1995).  The  input  variables,  called 
independent  variables,  are  under  the  control  of  the  researcher,  while  the  output  variables,  or 
the  response,  is  a  function  of  the  process  under  study.  This  concept  is  shown  graphically  in 

Figure  8,  where  |k  represents  the  process  inputs  and  E(y)  is  the  expected  value  of  the 

response,  y. 
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Process 


In  the  context  of  Figure  8  and  RSM,  the  MSS  model  is  the  process  under  study,  or  the 
black  box.  The  MSS  program  takes  in  a  number  of  variables  and  fits  either  a  diffusion 
coefficient  (Deff),  a  shape  factor  and  diffusion  path  length  (A.,  R),  or  gamma  function 
parameters  (a,  P)  by  minimizing  the  sum  of  squares  for  error  (SSE).  By  varying  the  input 
parameters,  this  research  examines  the  response  surface  of  SSE  with  the  original  radial 
geometry  function  and  with  a  gamma  distribution  for  fc  (8)  in  the  model  code  to  determine  if 
a  unique  solution  exists. 


!i  -► 

Inputs 

Process 

Output(s)  T|  =  E(y) 

Figure  8,  Graphical  Representation  of  Response  Process 


A  numerical  model  to  describe  sorption  and  desorption  with  perturbations  (liquid 
phase  added  or  removed  during  the  course  of  experimentation)  is  used  to  predict  sorption 
behavior  and  compare  predictions  of  the  MSS  model  to  the  data  sets.  The  governing 
equation  is  given  by  equation  (20) 

Vl®£+M^=0  (20) 

L  at  at 
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where  Vl  -  volume  of  liquid  phase, 

M  =  mass  of  paraffin,  and 

3S 

—  is  defined  by  the  MSS  model, 
dt 

At  some  time  t,  liquid  is  removed  and  may  be  added  to  the  reactor.  The  new 
boundary  conditions  are  defined  by 

Sj(t)  =  Sj.i(t) 

Vj -Cj (t)  =  Vj-i-Cj-1  (t)  -  VR(t)-Cj - 1  (t)  +  VA (t)-Cj (t) 

Vj-Vj.1  -  VR(t)  +  VA(t) 

where  j  =  boundary  condition, 

t  =  time  boundary  condition  changes, 

Vr  =  volume  of  liquid  removed, 

Va  =  volume  of  liquid  added,  and 

Cj(t)  =  concentration  of  liquid  during  boundary  condition  j  at  time  t. 

The  model  features  a  finite  difference  approach  employing  implicit  time  and  central 
space  differencing.  Two  different  versions  are  used,  featuring  the  composite  particle 
described  by  the  general  radial  geometry  function  and  by  the  gamma  density  function.  Both 
models  calculate  SSE  between  the  model  prediction  and  a  data  set,  and  are  included  in 
Appendix  C.  In  order  to  determine  the  nature  of  the  response  surface,  there  are  several 
phases  of  the  RSM  procedure  to  be  implemented. 

Screening 

Normally,  the  first  phase  of  RSM  involves  a  screening  process  to  determine  which  of 
the  independent  variables  significantly  influence  the  response.  In  this  research,  however, 
there  are  at  most  two  independent  variables  to  be  examined.  Therefore,  screening 
experiments  will  not  be  necessary.  The  radial  geometry  function  uses  X  and  R  as  the 
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independent  variables,  while  the  gamma  density  function  uses  a  and  p.  For  both  functions, 

the  response  is  SSE. 

Design  of  Experiments 

The  second  phase  of  RSM  uses  the  most  significant  variables  to  form  a  first-order 

approximation  of  the  SSE  response.  To  form  this  approximation,  factorial  designs  are 

typically  employed.  Factorial  designs  are  used  to  study  the  joint  effects  of  several 

independent  variables  on  a  response  (Montgomery,  1976).  By  using  two  levels  of  the 

independent  variables,  or  factors,  a  2k  factorial  design  is  created  where  there  are  k  total 

factors  being  studied.  The  factors  in  this  study,  X  and  R  or  a  and  p,  are  referred  to  as  factors 

2 

A  and  B,  respectively.  Each  factor  will  be  set  to  a  low  and  a  high  level,  resulting  in  a  2 
factorial  design,  where  all  possible  combinations  of  low  and  high  levels  result  in  four  runs,  or 
design  points.  Each  run  of  the  MSS  model  produces  a  response,  labeled  yi  through  y4,  as 
shown  in  Table  3. 

Table  3, 22  Factorial  Design 


Treatment 
Run#  Combination 

I 

Factorial  Effect 
A  B 

AB 

Response,  y 
SSE 

1  A  low,  B  low 

+ 

- 

- 

+ 

yi 

2  A  high,  B  low 

+ 

+ 

- 

- 

y2 

3  A  low,  B  high 

+ 

- 

+ 

- 

ys 

4  A  high,  B  high 

+ 

+ 

+ 

+ 

y4 

Using  the  designed  experiments  from  Table  3,  a  linear  approximation  of  the 
functional  relationship  between  A,  B,  and  y  is  produced.  This  is  accomplished  with 
regression  analysis,  a  statistical  technique  for  modeling  the  response  as  a  linear  combination 
of  various  forms  of  the  input  variables  and  their  interactions.  Before  performing  the 
regression  analysis,  the  input  variables  are  standardized,  or  coded,  using  equation  (21) 
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where  x ;  =  new  coded  variable, 

Xt  =  original  uncoded  variable, 

X  io  =  center  of  range  for  variable  X  ,• ,  and 
Si  =  half  of  the  range  for  variable  X 

This  transforms  the  variable  range  to  -1  and  +1,  corresponding  to  the  levels  under 
consideration,  low  and  high.  The  reason  for  coding  is  two-fold.  First,  it  increases  the 

2 

computational  ease  and  accuracy  of  the  model  coefficients  since  the  X’X  matrix  for  the  2 
design  will  have  some  of  the  off-diagonal  elements  equal  to  zero.  Second,  the  coefficient 
estimates  in  the  model  are  easier  to  interpret  because  the  parameter  estimates  of  the  linear 
model  are  uncorrelated  and  therefore  provide  statistically  independent  test  statistics  (Khuri 
and  Cornell,  1987). 

The  general  first-order  regression  model  is  y  =  X*P  +  £,  where  the  bold  symbols 
indicate  a  vector.  Errors  associated  with  this  model  fall  into  two  categories:  1)  random 
errors,  which  produce  the  variance,  a2,  arise  from  measurement  errors  of  the  process  under 
study  (included  in  e),  and  2)  systematic  errors,  or  bias  errors,  of  the  type  Xb  -  E(y)  (Myers 
and  Montgomery,  1995).  Since  the  MSS  model  is  deterministic,  there  are  no  random  errors, 
leaving  only  bias  errors  in  the  regression  model  (first-  or  second-order).  In  order  to  obtain 
the  minimum  bias,  a  minimum  bias  design  will  be  used  in  this  research. 

To  construct  a  minimum  bias  design,  the  factorial  design  points  are  pushed  closer  to 
the  center  of  the  design.  Using  the  coded  values  from  equation  (21),  the  new  points  would  be 
at  ±yj\p  instead  of  ±  1.  By  shrinking  the  overall  design  region,  model  bias  and  variance  are 

minimized  as  opposed  to  the  usual  criterion  of  minimizing  only  the  model  variance. 


Method  of  Steepest  Descent 

Using  the  first-order  response  model  of  SSE  obtained  through  regression  analysis,  a 
gradient  search  technique  called  the  method  of  steepest  descent  (MOSD)  iteratively 
determines  the  next  region  for  experimentation.  MOSD  is  continued  until  the  response 
(SSE)  no  longer  decreases  (or  decreases  only  slightly).  The  first  order  response  is  then 
augmented  with  interaction  effects  and  the  overall  model  fit  is  checked.  If  the  correlation 
coefficient  (R2)  is  not  satisfactory,  a  second  order  response  is  fit.  At  this  point  the  two-level 
factorial  design  is  augmented  to  include  center  points  in  order  to  test  for  curvature.  The 
central  composite  design  (CCD)  is  typically  used  for  fitting  a  second-order  response.  To 
make  the  CCD  a  minimum  bias  design,  the  axial  points  are  placed  at  ±  1  and  the  nonaxial 
points  are  placed  at  ±J\/2  .  Thus  all  design  points  are  maintained  on  a  circle  of  radius  equal 

to  one,  shrinking  the  design  region  under  investigation  and  reducing  the  bias. 

The  MSS  model  will  be  tested  with  the  radial  geometry  function  and  the  gamma 
density  function  to  get  a  picture  of  the  response  surface  in  the  minimum  SSE  region,  if  it 
exists.  The  hypothesis  is  that  changing  from  the  radial  geometry  function  to  the  gamma 
density  function  in  the  MSS  model  will  find  an  actual  minimum  SSE  rather  than  a  locus  of 
points  of  relatively  constant  SSE.  This  will  improve  the  capability  of  the  MSS  model  to 
predict  contaminant  concentrations  in  soil  systems.  A  notional  response  surface  is  given  in 
Figure  9,  which  shows  a  minimum  response  (at  least  locally). 
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IV.  RESULTS  AND  ANALYSIS 


Introduction 

The  purpose  of  this  thesis  is  to  improve  sorption  mass  transfer  models  by  using  a 
composite  geometry  to  describe  a  heterogeneous  mixture  of  sorbent  particles.  The 
improvement  comes  from  describing  a  composite  particle  geometry  using  a  statistical 
distribution  to  preserve  all  diffusion  path  lengths.  By  doing  so,  current  models  can  better 
predict  long-term  contaminant  residence  times. 

The  first  objective  is  to  attempt  to  explain  the  size  dependence  of  fitted  values  of  the 
diffusion  coefficient,  Deff,  from  paraffin  data  obtained  by  de  Venoge  (1996)  and  Mika 
(1997). 

The  second  objective  is  to  investigate  a  statistical  distribution  to  obtain  a  unique  fit  of 
two  parameters  that  describe  the  composite  geometry.  Mika  (1997)  found  that  fitting  the 
general  radial  geometry  function  yielded  nonunique  parameters.  The  measure  of  fit  is  sum  of 
squares  for  error  (SSE).  Can  a  minimum  SSE  be  found  for  a  pair  of  fitting  parameters,  or  is 
there  only  a  range  of  values  available? 

Size  Dependence  of  Derr 

The  MSS  model  predicts  different  diffusion  coefficients  in  two  separate  research 
efforts,  de  Venoge  (1996)  observed  that  the  fitted  Deff  depended  on  the  size  of  the  paraffin 
bead  being  modeled,  which  was  then  attributed  to  not  presoaking  the  paraffin  beads  in  a 
methanol/water  solution.  Mika  (1997)  presoaked  the  paraffin  beads  and  obtained  different 
diffusion  coefficients.  However,  these  new  coefficients  still  depend  on  the  paraffin  bead 
size,  de  Venoge  (1996)  shows  that  sphere  mass  is  normally  distributed  but  the  variance  of 
the  sphere  masses  is  dependent  on  the  sphere  size.  A  possible  explanation  is  that  the 
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variation  of  sphere  mass  (and  therefore  sphere  size)  causes  the  size  dependence  of  fitted  Deff 
coefficients. 

Mathematically,  the  quantity  D  a  /  5  describes  the  mass  transfer  rate  constant. 

Fitting  one  or  more  of  these  parameters  results  in  adjusting  parameters  so  D  a  /  5  remains 
constant.  For  a  given  solvent/sorbate/sorbent  system,  the  diffusion  coefficient  should  be 
constant  if  the  particle  size  and  shape  are  correctly  described. 

However,  if  the  geometry  (defined  by  the  size  index)  is  improperly  defined,  the  fitting 
program  adjusts  the  diffusion  coefficient  to  compensate  and  maintain  mass  balance.  In  this 
study,  beads  are  described  as  uniformly  sized  spheres  when  they  actually  should  have  been 
described  as  normally  distributed  sized  spheres. 

The  diffusion  coefficient  for  a  composite  shape  representing  normally  distributed 
spheres  can  be  estimated  by 


Dnd  =  D 


sphere 


(a/5) 

V  '  f  sphere 

Wl 


where  ND  =  normally  distributed. 

Using  the  sphere  statistics  in  Table  2,  a  distribution  of  radii  is  obtained  for  each 
sphere  size.  This  statistical  description  is  used  to  estimate  the  size  index  of  the  composite 
shape  of  the  large,  medium,  and  small  beads  using  the  procedure  for  a  homogeneous 
material,  heterogeneous  shape  described  in  Chapter  III. 

The  results  of  this  effort  are  presented  in  Table  4.  The  difference  in  size  indices 
between  the  composite  shape  and  the  nominal  uniform  sphere  shape  can  not  explain  the  size 
dependence  of  the  fitted  Deff  coefficients. 
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Table  4,  Diffusion  Coefficients  Adjusted  for  Bead  Size  Variance 

Bead  Fitted  Deff  Size  Index  Adjusted  Deff  Size  Index 
Size  (spherical)  (spherical)  (composite)  (composite) 


Large 

2.12 

136.0 

2.06 

139.8 

Medium 

2.01 

193.5 

1.94 

200.3 

Small 

1.96 

300.9 

1.89 

313.4 

Another  possible  explanation  for  the  size  dependence  of  Deff  has  to  do  with  the 
frequency  and  spacing  of  samples.  During  all  experiments,  samples  are  taken  at  the  same 
time  for  all  sizes  of  sorbents  within  a  given  condition.  But  since  the  smaller  beads  have  a 
faster  rate  of  sorption  (high  surface  to  volume  ratio),  the  samples  show  parts  of  the  rate  curve 
closer  to  equilibrium  than  would  the  same  samples  for  the  large  beads. 

In  order  to  examine  the  spacing  of  samples,  data  points  could  be  removed  from  a 
large  bead  data  set  before  using  it  in  the  MSS  model  to  calculate  Deff.  This  would  be  similar 
to  an  experiment  with  smaller  beads  where  samples  are  not  taken  until  the  system  is  closer  to 
equilibrium.  Figure  10  illustrates  the  effect  of  omitting  data  points  from  a  sorption 
experiment  by  Mika  (1997)  using  large  beads.  When  early  samples  are  deleted  from  the  data 
set,  the  fitted  Deff  coefficient  decreases.  Deletion  of  later  samples  has  little  effect  on  Deff 
indicating  the  importance  of  early  measurements  when  concentration  changes  rapidly  with 
time  and  where  most  of  the  overall  concentration  change  occurs.  This  is  the  factor  that 
explains  the  size  dependence  of  Deff.  It  is  therefore  important  to  capture  the  early 
concentration  changes  in  order  to  obtain  an  accurate  Deff. 
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Figure  10,  Loss  of  Data  Points  and  Relation  to  Deff 

RSM  and  Statistical  Distributions 

In  order  to  better  visualize  the  problem  of  fitting  nonunique  shape  and  size 
parameters  with  the  radial  geometry  function,  RSM  is  used  to  characterize  the  SSE  response 
surface  using  Mika’s  (1997)  even  size  bead  mix  data  set  found  in  Appendix  B.  A  detailed 
description  of  the  RSM  procedure  for  the  radial  geometry  function  is  in 
Appendix  D. 

Radial  Geometry  Function 

The  augmented  first-order  design  results  in  the  following  regression  equation  with 
R2  =  0.7945 

y  =  -  0.4195  +  0.2585  •  X  +  2.1298  •  R  -  1.21 65  •  X  •  R 
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Since  all  p-values  indicate  that  no  terms  are  significant,  a  second-order  response  surface  is 
fit.  The  resulting  regression  equation  is 

y  =  0.10856  +  0.05997*.  -  1.2371-R  +  0.0400  A2  -  0.9831  XR  +  6.71595R2 

with  a  correlation  coefficient  (R2)  of  0.9694.  Figure  1 1  is  a  plot  of  the  regression  equation. 
All  terms  are  left  in  the  equation  since  it  results  in  virtually  the  same  surface  as  with  only  the 
significant  terms.  Note  that  in  order  to  remain  in  the  trough  region  of  the  surface,  and 


Figure  11,  Radial  Geometry  Function  SSE  Response  Surface  From  Regression  Equation 

hence  near  a  minimum  SSE,  X  and  R  both  increase  as  the  trough  moves  away  from  the 
origin.  Figure  12  is  a  plot  of  the  true  surface  using  actual  data  from  the  MSS  model.  The 
figure  shows  the  nature  of  the  true  surface  encountered,  demonstrating  that  there  is  no  unique 


32 


best  fit  pair  of  X  and  R.  Both  figures  demonstrate  that  as  X  and  R  increase,  the  width  of  the 
trough  grows.  This  is  another  limitation  of  the  radial  geometry  function,  since  real  soil 
systems  are  probably  characterized  by  composite  shape  factors,  X  »  2.  Figure  13  is  a  plot  of 
log(SSE)  contours,  allowing  a  better  visualization  of  the  trough.  Figure  13  shows  that  the 
trough  widens  as  the  surface  extends  from  the  origin.  Although  not  visible  from  this  plot,  the 
SSE  values  also  decrease  slightly  as  the  trough  widens. 


l.°°1 


Figure  12,  Radial  Geometry  Function  SSE  Surface  From  MSS  Model  Data 


X 


Figure  13,  Radial  Geometry  Function  Contour  Plot  of  log  (SSE) 
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Gamma  Density  Function 

After  changing  the  radial  geometry  function  in  the  MSS  model  code  to  a  gamma 
density  function,  Mika’s  (1997)  even  size  bead  mix  data  set  is  again  used  to  evaluate  the  SSE 
response  with  RSM.  By  using  a  gamma  density  function  instead  of  the  radial  geometry 
function,  a  unique  pair  of  fitting  parameters  is  found.  To  correspond  with  the  gamma 
function,  these  parameters  are  named  a  and  p.  The  shape  parameter,  a,  controls  the  shape  of 
the  gamma  function.  Increasing  it  has  an  effect  similar  to  decreasing  X  in  the  radial  geometry 
function.  The  scale  parameter,  p,  is  similar  to  R  in  the  radial  geometry  function,  although  P 
is  dimensionless.  The  radial  geometry  function  is  intuitively  appealing  since  R  can  be 
visualized  as  the  actual  particle  radius  in  centimeters  and  X  =  0,  1 ,  or  2  represents  a  layer, 
cylinder,  or  sphere,  respectively,  as  mentioned  in  Chapter  II.  Unfortunately,  a  and  P  do  not 
have  this  same  physical  interpretation. 

The  first-order  augmented  design  for  the  gamma  function  MSS  model  is 
y = 0.7476  -  0.8629  •  a  -  8.4792  •  P  +  9.860  •  a  •  P 
with  a  corresponding  R2  of  0.9271 .  At  first  glance,  this  model  appears  to  be  adequate. 
However,  the  p-values  for  all  parameter  estimates  are  >  0.293,  indicating  a  better  model  is 
available.  The  second-order  response  function  is 

y  =  1.5096  - 1.9220  •  a  -  15.6135  •  P  +  0.6155  •  a2  +  9.8676  •  a  •  P  +  41.4526  •  p2 
where  R2  =  0.9999  and  all  p-values  are  <  0.0001.  Figure  14  is  a  plot  of  the  SSE  values 
calculated  by  the  MSS  model  using  the  gamma  density  function.  A  plot  of  the  above 
regression  equation  looks  similar,  but  does  not  show  the  minimum  as  clearly  as  the  plotted 
values  of  a,  P,  and  SSE  in  Figure  14.  This  plot  clearly  demonstrates  that  a  unique  minimum 
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pair  of  fitting  parameters  exists  for  the  MSS  model  when  using  the  gamma  density  function. 
The  area  of  the  minimum  seems  quite  large  in  Figure  14  since  0.84  <  a  <  0.95  and  0.075  <  3 
<  0.095  will  result  in  approximately  the  same  SSE.  However,  a  and  P  can  range  from  zero  to 
infinity  for  the  gamma  function.  Therefore,  the  area  of  the  minimum  in  Figure  14  is 
relatively  small. 


Figure  1 4,  Gamma  Density  Function  SSE  Surface  From  MSS  Model  Data 

There  is  a  definite  advantage  of  the  gamma  MSS  model  over  the  radial  geometry 
version  of  the  MSS  model.  The  gamma  MSS  model  will  obtain  a  unique  pair  of  fitting 
parameters  for  a  given  soil  mixture.  The  radial  geometry  MSS  model  will  not  fit  unique 
pairs  of  parameters,  only  a  unique  size  index.  If  unique  parameters  for  both  size  and  shape  of 
composite  sorbent  particles  can  be  fit  from  data,  these  parameters  can  characterize  natural 
sorbents  (soil  systems)  more  thoroughly  than  the  single  parameter,  size  index.  Even  if  the 
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size  and  shape  parameter  have  no  physical  interpretation  as  in  the  radial  geometry  function, 
different  soil  systems  can  be  compared  based  on  unique  fitting  parameters. 

Real  soil  systems  are  heterogeneous  mixes  of  particle  sizes  that  are  only  roughly 
approximated  by  even  or  uneven  mixes  of  paraffin  beads.  Since  the  even  mix  of  paraffin 
beads  in  Mika’s  (1997)  data  has  a  unique  minimum  pair  of  fitting  parameters,  the  obvious 
test  of  the  new  gamma  MSS  model  is  to  see  if  it  will  find  a  unique  pair  of  fitting  parameters 
for  a  real  soil  system. 

Borden  Data  Set 

The  Borden  data  set  as  described  in  Chapter  III  is  based  on  real  soil.  Appendix  A 
contains  a  simulated  data  set  for  the  MSS  model  based  on  the  Borden  soil.  Using  this 
simulated  data,  the  RSM  procedure  is  applied  to  the  gamma  function  version  of  the  MSS 
model  to  see  if  it  will  find  a  unique  pair  of  fitting  parameters  for  real  soil. 

The  augmented  first-order  design  for  the  Borden  data  is 

y  -  0.561 8  -  3.75 1 2  •  a  -  0.1 993  •  P  +  3.0773  •  a  •  P 
The  correlation  coefficient  for  this  equation  is  only  0.007,  while  the  unaugmented  first-order 
design  has  R2  =  0.8508.  The  lower  R2  is  due  to  the  curvature  of  the  response  surface  and 
indicates  the  need  for  a  second-order  response  to  be  fit.  The  second-order  response  surface 
equation  is 

y  =  - 13.193- 12.1518-a  + 16.854  P  +  44.7592-a2  +2.7823-a-P-5.2112-p2 
with  R2  =  0.9137.  Figure  15  is  a  plot  of  the  MSS  model  SSE  responses  with  the  Borden  data 
set. 

For  a  real  soil  data  set,  the  gamma  version  of  the  MSS  model  finds  an  actual 
minimum  pair  of  fitting  parameters.  From  Figure  15  it  appears  that  there  are  actually  two 
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regions  where  a  minimum  can  occur.  However,  the  actual  data  points  show  that  the  true 
minimum  is  within  the  larger  of  the  two  circular  areas.  The  figure  just  doesn’t  include 
enough  lines  to  show  this  fact.  The  significant  point  is  that  the  gamma  MSS  model  locates  a 
unique  minimum  for  real  soil. 


Figure  15,  Borden  Soil  SSE  Surface  From  MSS  Model  Data 


With  the  Borden  soil  data  set,  a  «  0.1  in  the  gamma  MSS  model,  and  X  «  8.0  in  the 
radial  geometry  MSS  model.  This  value  of  a  is  very  low  considering  that  there  is  no  upper 
limit  to  the  range  for  a.  For  Mika’s  (1997)  even  mix,  a  «  0.9  in  the  gamma  MSS  model  and 
the  corresponding  shape  factor  in  the  radial  geometry  function  is  approximately  2.7.  The 
significance  is  that  in  the  radial  geometry  function,  the  shape  factor  increases  for  real  soil  and 
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has  the  range  to  do  so.  For  the  gamma  MSS  model,  this  range  is  not  available  since  the 
shape  factor  decreases  for  real  soil  and  has  a  lower  limit  of  zero.  The  small  available  range 
for  the  shape  factor  in  the  gamma  MSS  model  may  be  the  reason  that  this  model  is  able  to 
locate  unique  fitting  parameters.  Further  investigation  is  necessary  in  order  to  determine 
whether  or  not  this  is  true. 

A  potential  problem  with  the  gamma  MSS  model  is  that  low  a  values  are  difficult  to 
use  in  the  gamma  function.  The  gamma  MSS  model  is  only  calibrated  for  a  >  0.35.  Before 
actually  characterizing  soils  with  the  gamma  MSS  model,  more  sensitivity  analysis  and 
model  calibration  are  required. 

Since  the  gamma  function  is  mathematically  difficult  to  use,  the  lognormal  or 
Weibull  density  functions  might  be  used  instead.  The  lognormal  function  has  been  used  by 
Heyse  (1994)  and  Culver  et  al.  (1997)  to  describe  mass  transfer  rates.  Both  the  lognormal 
(equation  10)  and  Weibull  (equation  11)  density  functions  are  two-parameter  distributions 
that  can  assume  a  wide  variety  of  shapes  based  on  input  parameters.  These  distributions 
have  the  potential  to  be  as  useful  as  the  gamma  distribution  described  in  this  thesis. 
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V.  SUMMARY  AND  CONCLUSIONS 

Two  problems  found  in  earlier  research  with  the  MSS  model  are  examined  in  this 
thesis.  Both  de  Venoge’s  (1996)  and  Mika’s  (1997)  research  find  that  fitted  Deff  coefficients 
depend  on  sphere  sizes.  The  variations  in  Deff  are  not  due  to  the  normally  distributed  sphere 
masses  (and  hence  sphere  sizes),  but  to  the  spacing  of  sampling  times  during  sorption  / 
desorption  experiments.  Most  of  the  concentration  changes  occur  early  in  the  sorption 
experiment.  These  changes  must  be  captured  to  accurately  describe  the  sorption  rate  curve 
and  fit  model  parameters. 

When  using  the  MSS  model  to  fit  composite  particle  geometry  parameters,  nonunique 
fits  are  obtained.  This  is  due  to  limitations  of  the  radial  geometry  function  (equation  8)  as 
the  shape  factor  increases.  When  the  gamma  density  function  (equation  9)  is  coded  into  the 
MSS  model,  unique  fits  are  obtained  for  synthetic  and  natural  sorbent  data  sets. 

In  the  course  of  determining  the  answers  to  these  problems,  several  factors  became 
apparent  that  merit  further  attention.  First,  to  improve  future  sorption  experiments,  more 
frequent  sampling  is  needed  as  the  bead  size  decreases.  This  will  capture  earlier  portions  of 
the  sorption  process  previously  lost  due  to  the  faster  approach  to  equilibrium  of  the  high 
surface  area  to  volume  ratio  spheres. 

The  radial  geometry  MSS  model  is  an  improvement  over  the  equivalent  sphere  model 
since  it  preserves  the  longer  path  lengths.  The  gamma  MSS  model  is  an  even  better 
improvement  since  it  will  locate  unique  best  fit  geometry  parameters  for  composite  sorbents, 
whether  that  sorbent  is  natural  soil  or  synthetic  material.  Continued  investigation  is  needed, 
however,  to  calibrate  the  model  for  a  <  0.35  since  the  shape  parameter  takes  on  values 
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substantially  below  0.35  for  natural  soils.  A  better  algorithm  for  computing  the  gamma 
function  may  simplify  the  calibration  effort. 

Finally,  since  the  gamma  function  is  difficult  to  evaluate,  the  lognormal  and  Weibull 
density  functions  should  be  examined.  These  functions  offer  a  potential  advantage  over  the 
gamma  function  because  they  are  computationally  more  efficient  than  the  gamma  density 
function.  If  these  functions  can  also  fit  unique  geometry  parameters,  characterization  of  soil 
types  can  be  accomplished. 


APPENDIX  A:  Borden  Soil  Data  Set 


Borden  Soil  -  Gamma  MSS  model  simulation:  Run  5 
RKD  DEFF  DPATH  RLAM  XSSO  FAST 

0.727  0.3390E-10  0.011  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 
0  1  10  0  0 
1 

0.013481889  0.09779271 
CONST ANTS(Ms,  VLi,  CLi) 

1 

0.10000E+01  0.10000E+01  0.10000E+01  0.10000E+01  10 
1  0.10000E+03  0.00000E+00  0.00000E+00  0.00000E+00 

1  0.10000E+04  0.00000E+00  0.00000E+00  O.OOOOOE+OO 

1  0.10000E+05  O.OOOOOE+OO  0.00000E+00  0.00000E+00 

1  0.10000E+06  0.00000E+00  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.10000E+07  O.OOOOOE+OO  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.20000E+07  O.OOOOOE+OO  0.000Q0E+00  O.OOOOOE+OO 

1  0.30000E+07  O.OOOOOE+OO  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.40000E+07  O.OOOOOE+OO  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.50000E+07  O.OOOOOE+OO  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.60000E+07  O.OOOOOE+OO  O.OOOOOE+OO  O.OOOOOE+OO 

1  0.001 

EXP  TIME(SEC)  CONCENTRATION(C) 

1  0.640000000E+03  0.967326343E+00 

1  0. 1 90000000E+04  0.944389999E+00 

1  0.640000000E+04  0.903672993E+00 

1  0. 1 90000000E+05  0.84875 1038E+00 

1  0.550000000E+05  0.778037339E+00 

1  0.9 1 0000000E+05  0.741 21 8477E+00 

1  0. 1 90000000E+06  0.689654505E+00 

1  0.280000000E+06  0.66561 8 181E+00 

1  0.460000000E+06  0.638607705E+00 

1  0.640000000E+06  0.623 186839E+00 

1  0. 820000000E+06  0.613297594E+00 

1  0. 1 1 0000000E+07  0.603763425E+00 

1  0. 1 30000000E+07  0.599406350E+00 

1  0. 1 80000000E+07  0.59282541 3E+00 

1  0.220000000E+07  0.589753509E+00 

1  0.260000000E+07  0.587632942E+00 

1  0.350000000E+07  0.584662259E+00 

1  0.440000000E+07  0.582938743E+00 

1  0.520000000E+07  0.581937766E+00 

1  0.6 1 0000000E+07  0.581267679E+00 

1  0.690000000E+07  0.580767 190E+00 
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APPENDIX  B:  Mika’s  (1997)  Even  Data  Set,  Condition  5 

mharrison  TEST  FOR  MSS  model,  Run  1 :  even,  two  conditions 
RKD  DEFF  DPATH  RLAM  XSSO  FAST 

9.35  2.0301  E-08  0.2453  2.0  0.0  0.0 

1  MEANS  FIT  PARAMETER,  0  MEANS  KEEP  IT  FIXED 

0  0  1  0  0  0 

9 

1  0.1 

2  0.1 

1  0.25 

2  0.25 

1.5  0.175 

1.574071  0.249172 

1.648142  0.323345 

1.722213  0.397517 

1.796284  0.471690 

CONST ANTS(Ms,  VLi,  CLi) 

2 

2.02006  20.33937  0.979  1.0  28 

2.02649  20.40128  0.977  1.0  28 


1 

69840.00 

0.56231 

0. 

0. 

1 

166560.0 

0.45357 

0. 

0. 

1 

251460.0 

0.41852 

0. 

0. 

1 

413160.0 

18.38716 

20.11721 

0. 

1 

585060.0 

0.36828 

0. 

0. 

1 

753960.0 

0.36446 

0. 

0. 

1 

948360.0 

0.41360 

0. 

0. 

1 

1277760. 

19.07225 

20.11115 

0.997 

1 

1361880. 

0.41300 

0. 

0. 

1 

1468560. 

0.40200 

0. 

0. 

1 

1549860. 

0.43592 

0. 

0. 

1 

1632840. 

0.45680 

0. 

0. 

1 

1793760. 

0.37290 

0. 

0. 

1 

2921640. 

18.01844 

20.02739 

0. 

1 

3022440. 

0.39692 

0. 

0. 

1 

3102840. 

0.40276 

0. 

0. 

1 

3193440. 

0.38009 

0. 

0. 

1 

3275040. 

0.40635 

0. 

0. 

1 

3371040. 

0.46096 

0. 

0. 

1 

3448560. 

0.44722 

0. 

0. 

1 

3610320. 

17.57318 

20.06768 

0.997 

1 

3699960. 

0.42443 

0. 

0. 

1 

3786060. 

0.39373 

0. 

0. 

1 

3873060. 

0.46827 

0. 

0. 

1 

3963060. 

0.46509 

0. 

0. 
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1 

4043460. 

0.43136 

0. 

0. 

1 

4134660. 

0.41749 

0. 

0. 

1 

4306260. 

0.35215 

0. 

0. 

2 

69840.00 

0.43149 

0. 

0. 

2 

166560.0 

0.35314 

0. 

0. 

2 

251460.0 

0.35427 

0. 

0. 

2 

413160.0 

18.86931 

19.96301 

0. 

2 

585060.0 

0.36885 

0. 

0. 

2 

753960.0 

0.39213 

0. 

0. 

2 

948360.0 

0.39560 

0. 

0. 

2 

1277760. 

18.74505 

20.28702 

0.997 

2 

1361880. 

0.38462 

0. 

0. 

2 

1468560. 

0.38112 

0. 

0. 

2 

1549860. 

0.48235 

0. 

0. 

2 

1632840. 

0.42297 

0. 

0. 

2 

1793760. 

0.29567 

0. 

0. 

2 

2921640. 

18.45080 

19.91980 

0. 

2 

3022440. 

0.38641 

0. 

0. 

2 

3102840. 

0.42475 

0. 

0. 

2 

3193440. 

0.39736 

0. 

0. 

2 

3275040. 

0.40122 

0. 

0. 

2 

3371040. 

0.41289 

0. 

0. 

2 

3448560. 

0.48550 

0. 

0. 

2 

3610320. 

17.45752 

20.07969 

0.997 

2 

3699960. 

0.43547 

0. 

0. 

2 

3786060. 

0.40964 

0. 

0. 

2 

3873060. 

0.42862 

0. 

0. 

2 

3963060. 

0.42927 

0. 

0. 

2 

4043460. 

0.36309 

0. 

0. 

2 

4134660. 

0.42079 

0. 

0. 

2 

4306260. 

0.43494 

0. 

0. 

1  0.001 

EXP 

TIME(SEC) 

CONCENTRATION(C) 

1 

166560.0 

0.599 

1 

251460.0 

0.570 

1 

413159.0 

0.540 

1 

413161.0 

0.014 

1 

585060.0 

0.209 

1 

753960.0 

0.227 

1 

948360.0 

0.234 

1 

1277759. 

0.239 

1 

1277761. 

0.982 

1 

1361880. 

0.739 

1 

1468560. 

0.676 

1 

1549860. 

0.661 

1 

1632840. 

0.649 

1 

1793760. 

0.641 

1 

2921639. 

0.624 

1 

2921641. 

0.013 

1 

3022440. 

0.201 

1 

3102840. 

0.230 

1 

3193440. 

0.257 

1 

3275040. 

0.280 

1 

3371040. 

0.291 

1 

3448560. 

0.293 

1 

3610319. 

0.302 

1 

3610321. 

0.984 

1 

3699960. 

0.765 

1 

3786060. 

0.721 

1 

3873060. 

0.705 

1 

3963060. 

0.691 

1 

4043460. 

0.686 

1 

4134660. 

0.680 

1 

4306260. 

0.667 

1 

5085960. 

0.651 

2 

69840.00 

0.675 

2 

166560.0 

0.599 

2 

251460.0 

0.571 

2 

413159.0 

0.538 

2 

413161.0 

0.010 

2 

585060.0 

0.211 

2 

753960.0 

0.229 

2 

948360.0 

0.235 

2 

1277759. 

0.241 

2 

1277761. 

0.981 

2 

1361880. 

0.740 

2 

1468560. 

0.686 

2 

1549860. 

0.663 

2 

1632840. 

0.650 

2 

1793760. 

0.642 

2 

2921639. 

0.626 

2 

2921641. 

0.010 

2 

3022440. 

0.202 

2 

3102840. 

0.230 

2 

3193440. 

0.259 

2 

3275040. 

0.280 

2 

3371040. 

0.291 

2 

3448560. 

0.290 

2 

3610319. 

0.301 

2 

3610321. 

0.988 

2 

3699960. 

0.767 

2 

3786060. 

0.724 
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2  3873060.  0.705 

2  3963060.  0.690 

2  4043460.  0.689 

2  4134660.  0.679 

2  4306260.  0.672 

2  5085960.  0.653 
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APPENDIX  C:  MSS  Models 


Radial  Geometry  Function  Version 

C 

C  Program  Testfit 

C 

C  Adapted  for  batch  system  experiencing: 

C  Multisite  (series)  sorption  on  soil  (variable  distribution) 

C  Liquid-phase  added/removed  perturbations 

C  Simultaneous  Fitting  multiple  experiments 

C 

C  THIS  VERSION  CONFIGURED  TO  CALCULATE 

C  SSQ  FOR  ERROR  ONLY 

C 

C  NUMERICAL  SOLUTION 

C 

C  BY:  E.  Heyse  and  D.  C.  Coulliette 
C  University  of  Florida,  13  Feb  97 

C 

C  dim  scrat  =  5*np  +np*np  +(2+np)*nob 

C  signs  =  1  means  parms  cannot  change  sign 
C  P(IPOINT(I))  =  VAL(I)  —  PNAME(I) 

C 

C  IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  RKD,DEFF,DPATH,RLAM,XSSO,FAST 
DOUBLE  PRECISION  XM(10),V(10),PCO(10),SFO(10),XEFT(10),TF 
DOUBLE  PRECISION  TBCC(220),VR(220),VA(220),CA(220) 
DOUBLE  PRECISION  X(250) 

DOUBLE  PRECISION  ANS(0:500001),TIME,CTIME,STOLD,RSLOW 
DOUBLE  PRECISION  XIC(501),DELTA(501),FN(501),CL,VS 
DOUBLE  PRECISION  XLAM(50),XR(50),SSQ 
INTEGER  FDOMAIN 
DIMENSION  IX(250),Y(250) 

DIMENSION  IFC(220)  ,NEFD(  1 0) 

CHARACTER*  8  PNAME(7) 

CHARACTER*  15  INFIL,OUTl 
DIMENSION  V  AL(7)  ,IPOINT  (7)  ,1 V  ARY  (7) 

DIMENSION  F(250) 

COMMON  /BMOD/XM,V,PCO,SFO,XEFT,NEFD,TBCC,VR,VA,CA 
COMMON  /BMODA/IFLAG,TF 

COMMON  /VALUE/  VAL,IPOINT,IVARY,PNAME,NFIX 
COMMON  /BDAT/  IX,X,NEXP 
CHARACTER*  80  TIT, LIN  1  ,LIN2,LIN3,LIN4 

nfix=6 

PNAME(1)='  Kd ' 


PNAME(2)='  Deff ' 

PNAME(3)='  dpath ' 

PNAME(4)=' Lambda' 

PNAME(5)='  SsO ' 

PNAME(6)='  Fast ' 

WRITE(*,*)  'NAME  OF  INPUT  FILE?' 

READ(5,1 1 1 1)INFIL 
1111  FORMAT(A15) 

WRITE(6,*)'NAME  OF  OUTPUT  FILE?' 

READ(5,1 1 1  l)OUTl 

OPEN(UNIT=3  ,FILE=INFIL,ST  ATUS='OLD') 

OPEN(UNIT=l  0,FILE=OUT1  ,STATUS='NEW') 

READ(3,'(A)')  TIT 
C 

C  Read  parameters  to  be  fitted 
C 

C  RKD  =  partition  coefficient,  KP, 

C  DEFF  =  Coefficient  of  Diffusion,  D_eff  cmA2/ sec 

C  DPATH  =  Total  diff.  path  length,  d,  cm 
C  RLAM  =  shape  factor  Lambda 

C  XSSO  -  Initial  Concentration  in  slow  soil  sites,  Sso,  (Xg/g 
C  FAST  =  Fraction  of  equilibrium  sorption  sites 
C 

READ(3,'(A)')  LIN1 

READ(3,*)  RKD, DEFF, DPATH, RLAM, XSSO, FAST 
C 

C  Read  Constants 

C 

C  FOC  =  Mass  fraction  organic  carbon  on  soil,  foe 
C  XM  =  Mass  of  soil  in  reactor,  g 
C  V  =  Volume  of  solvent  in  reactor  at  time  zero,  mL 
C  PCO  =  initial  concentration  in  liquid,  mg/1 

C  SFO  =  experiment  specific  initial  concentration  in  solid 
C 

C  Note:  Intial  cone  in  solid  phase  =  XSS0*SF0(IEXP) 

C  if  initial  sorbed  cone  is  the  same  for  all  experiments, 

C  make  XSSO  =  S(t=0)  and  SFO(IEXP)  =  1 .0 

C  if  initial  sorbed  cone  is  different  for  each  experiment, 

C  make  XSSO  =  1 .0  and  SFO(IEXP)  =  S(t=0,  exp=IEXP) 

C 

C  NEFD  =  Number  of  boundary  condition  changes 
C  IFC  =  boundary  condition  number 
C  TBCC  =  time  for  boundary  condition  to  end 
C  XQ  =  Flow  rate,  mL/min 
C  XCIN  =  Dimensionless  inlet  concentration 
C  XMR  =  Mass  of  reactor,  g 
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C  XKPR  =  Linear  partition  coef.  on  reactor,  mL/g 

C  XKR  =  1  st  order  mass  transfer  rate  onto  reactor,  minA- 1 

C  XFR  =  Fraction  instantaneous  sorbing  reactor  sites 
C  VS(J)  =  volume  of  liquid  in  rx  for  BC  period  j,  ml 

C  VR(J)  =  volume  of  liquid  removed  at  end  of  BC  period  J,  ml. 

C  VA(j)  =  volume  of  liquid  added  at  end  of  BC  period  j,  ml. 

C  CA(J)  =  concentration  in  liquid  added  at  end  of  BC  period  J,  mg/1. 

C  XIC  =  sorbed  concentration  in  each  slow  compartment,  |ig/g 
C  CL  =  current  liquid  concentration,  mg/1 
C 

READ(3,'(A)')  LIN2 
READ(3,*)  (IVARY(I),I=1  ,NFIX) 

READ  IN  HIGH  AND  LOW  values  for  lambda  and  R 


READ(3,*)NPAIR 

DO  7  IPAIR=1,NPAIR 

READ(3,*)  XLAM(IPAIR),  XR(IPAIR) 

7  CONTINUE 

READ(3,'(A)')  LIN3 
READ(3,*)NEXP 
NABC=0 
DO  17  I=1,NEXP 

READ(3,*)  XM(I),V (I),PCO(I),SFO(I),NEFD(I) 

C  READ(3,*)  XMR, XKPR, XKR, XFR 
NABC=NABC+NEFD(I) 

17  CONTINUE 

DO  6  I=1,NABC 

6  READ(3,*)  IFC(I) ,TBCC(I) , VR(I) , V A(I) ,C A(I) 

C  Read  numerical  parameter 

C  IFLAG  =  0  for  sorbed/liquid  concentration 
C  IFLAG  =  1  for  liquid  concentration 
C  TF  =  time  factor  =  time  of  time  step/time  for  mass  transfer 
C  Try  TF  =  0.01 

READ(3,*)  IFLAG,  TF 
C 

C  Read  Data  to  be  fitted 

C 

C  X  =  time,  seconds 

C  Y  =  observation  at  time  X,  see  IFLAG. 

C 

C  XEFT(I)  =  experiment  final  time  (last  observation  for  each  exp.) 

C 

C 

READ(3/(A)')  LIN4 
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NOB=0 
DO  3  1=1,250 

READ(3,*,END=4)  IX(I),X(I),Y(I) 

IF  ((I.GE.2).AND.(IX(I).NE.IX(I- 1 )))  XEFT(IX(I-1))=X(I-1) 

3  NOB=NOB+l 

4  CONTINUE 
XEFT(NEXP)=X(NOB) 

C 

C  Batch  model  for  use  with  fitting  routine  TESTFIT 

C 

C  Set  up  to  handle  multiple  experiments  (10) 

C 

C  AT  THIS  POINT  WE  CYCLE  THROUGH  VARIOUS  VALUES  OF 

C  PARAMETERS  TO  CACULATE  SSQ  BETWEEN  OBSERVED 

C  AND  PREDICTED  RESPONSES 

C 

DO  300  IPAIR2=  1  ,NPAIR 
RLAM=XLAM(IPAIR2) 

DPATH=XR(IPAIR2) 

C 

VAL(1)=RKD 
VAL(2)=DEFF 
V  AL(3)=DPATH 
VAL(4)=RLAM 
VAL(5)=XSS0 
VAL(6)=FAST 
NFRAC=(RLAM+2.dO)*  10+1 
IF  (NFRAC.LT.50)  NFRAC=50 
IF  (NFRAC.GT.500)  WRITE(6,*)  'NFRAC  = NFRAC 
IF  (NFRAC.GT.500)  WRITE(10,*)  'NFRAC  =  NFRAC 
IF  (NFRAC.GT.500)  WRITE(4,*)  'NFRAC  = ',  NFRAC 


Loop  through  each  experiment  and  set  up  appropriate 
initial  conditions 

NSUM=0 

IT=1 

DO  200  IEXP=1,NEXP 

CL=(PC0(IEXP)  *  V  (IEXP)+XM(IEXP)  *  XSS0*  SFO(IEXP)  *FAST)/ 
.  (V(IEXP)+RKD*XM(IEXP)*FAST) 

VS=V(IEXP) 


Input  the  domain  partitioning.  Delta(I)  is  the  delta  for 
each  compartment. 
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DO  10  I=1,NFRAC 
DELTA(I)-DPATH/NFRAC 
10  CONTINUE 
C 

C  Subroutine  INTF  generates  the  F(i)'s  by  integrating  f(delta) 
C  with  RLAM  as  the  shape  parameter 
C 

C  Set  the  initial  condition  within  the  immobile  zone 
C 

CALL  INTF(RLAM,  DELTA,  DPATH,NFRAC,FN) 

C 

C  STOLD  is  the  sum  of  the  S  for  the  diffusion  compartments 
C 

STOLD=0.0 
DO  20  I=1,NFRAC 

XIC(I)=XSS0*SF0(IEXP)*FN(I)*(1  .DO-FAST) 
STOLD=STOLD+XIC(I) 

20  CONTINUE 
C 

C  Loop  through  each  flow  domain 
C 

FDOMAIN=l 

TIME=0.0 

C  Note  the  final  time  is  X(NOB) 

DO  WHILE  (FDOMAIN.LE.(NEFD(IEXP)-l- 1 )) 

C 

C  Calculate  ANS  for  beginning  of  boundary  condition 

C 

IF  (IFLAG.EQ.l)  THEN 
ANS(0)=CL 

ELSE 

ANS(0)=(STOLD+FAST*RKD*CL)/CL 
END  IF 
C 

C  CTIME  is  the  cumulative  solution  time 
CTIME=TBCC(FDOMAIN+NSUM) 

IF  (FDOMAIN.EQ.(NEFD(IEXP)+l ))  THEN 
DO  210  JCOUNT=l,NOB 
IF  (IX(JCOUNT).EQ.IEXP)  CTIME=X(JCOUNT) 

210  CONTINUE 
END  IF 
C 

C  Calculate  retardation  factor  for  slow  sites 
C 

RSLOW=(XM(IEXP)*RKD*(l.dO-FAST))/VS 

C 


50 


C  Calculate  the  number  of  time  steps  for  the  BC 

C 

ITHALF=(CTIME-TIME)*RSLOW*FN(  l)*DEFF/(TF*DELTA(l)**2.dO)+l 
ITSTEPS=ITHALF*2 
C 

C  Set  minimum  value  of  ITSTEPS 
C 

IF  (ITSTEPS.LT.100)  GOTO  30 
GOTO  35 

30  IF  ((CTIME-TIME).GT.(XEFT(IEXP)/200.d0))  ITSTEPS=100 

C  IF  ((ITSTEPS.LT.  1 00).  AND.((CTIME-TIME).GT.(XEFT(IEXP)/ 

C  200.d0)))  ITSTEPS=100 

C 

C  Warn  if  number  of  time  steps  exceeds  array  size 

C 

35  CONTINUE 

IF  (ITSTEPS.GT.500000)  WRITE(6,*)  'ITSTEPS  = ITSTEPS 
IF  (ITSTEPS.GT.500000)  WRITE(10,*)  'ITSTEPS  = ',  ITSTEPS 
C  IF  (ITSTEPS.GT.500000)  WRITE(4 ,*)  'ITSTEPS  - ',  ITSTEPS 

IF  (ITSTEPS.GT.500000)  GOTO  300 
C 

C  FDSOLVE  computes  the  soln  from  time  to  crime  for  itsteps  and 
C  passes  the  soln  back  as  array  ANS(ITSTEPS) 

C 

CALL  FDSOLVE(TIME,CTIME,  ITSTEPS,  ANS,XIC, 
FN,DELTA,CL,STOLD,VS,IEXP,NFRAC) 

C 

C  INTERP  interpolates  the  soln  from  FDSOLVE  onto  the  observation 
C  times  and  stores  the  results  in  array  F(nobs) 

C 

DO  WHILE  ((X(IT).LE.CTIME).AND.(X(IT).GE.TIME) 

.AND.(IX(IT).EQ.IEXP)) 

CALLINTERP(X(IT),TIME,CTIME,ANS,ITSTEPS,F(IT)) 

IT=IT+1 
END  DO 
TIME=CTIME 
C 

C  Adjust  aqueous  volume  and  new  aqueous  cone. 

C 

RFASTV =XM(IEXP)*RKD*FAST +V  S 
XMADD=VA(FDOMAIN+NSUM)*CA(FDOMAIN+NSUM) 
XMREM=VR(FDOM  AIN+NSUM)  *CL 
XNEWM=RFASTV*CL+XMADD-XMREM 

CL=XNEWM/(RFASTV-VR(FDOMAIN+NSUM)+VA(FDOMAIN+NSUM)) 

C  CL=(VS*CL+VA(FDOMAIN+NSUM)*CA(FDOMAIN+NSUM)- 

C  VR(FDOM  AIN +N SUM)  *  CL)/ 
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C  (VS-VR(FDOMAIN+NSUM)+VA(FDOMAIN+NSUM)) 

VS=VS-VR(FDOMAIN+NSUM)+VA(FDOMAIN+NSUM) 

FD0MAIN=FD0MAIN+1 

C 

END  DO 

NSUM=NSUM+NEFD(IEXP) 

200  CONTINUE 
C 

SSQ=0.d0 
DO  253  IS=l,NOB 
SSQ=SSQ+(Y(IS)-F(IS))**2.dO 
253  CONTINUE 

WRITE(10,*)  RLAM,DPATH,SSQ 
300  CONTINUE 
C  Output  the  results 
C 

CDLC  1=1 

CDLC  DO  WHILE  (I.LT.NOB) 

CDLC  WRITE(*,*)I,F(I) 

CDLC  1=1+1 

CDLC  END  DO 

CDLC  STOP 
2222  CONTINUE 
STOP 
END 
C 

SUBROUTINE  FDSOLVE(TIME,CTIME,ITSTEPS,ANS,XIC, 

.FN,  DELTA,  CL,  STOLD,VS,IEXP,NFRAC) 

C 

C  IMPLICIT  NONE 

DOUBLE  PRECISION  RKD,DEFF,DPATH,RLAM,XSSO,FAST 
DOUBLE  PRECISION  XM(10),V(10),PC0(10),SF0(10),XEFT(10),TF 
DOUBLE  PRECISION  TBCC(220),VR(220),VA(220),CA(220) 
DOUBLE  PRECISION  ANS(0: 500001), TIME, CTIME,STOLD 
DOUBLE  PRECISION  XIC(501),DELTA(501),FN(501),CL,VS 
DOUBLE  PRECISION  K(501),RES(501),A(501),B(501),RESN,ANS2 
INTEGER  ITSTEPS  ,NFRAC 

DOUBLE  PRECISION  SNEW(501),SOLD(501),ST,DELTAT 
DIMENSION  IFC(220)  ,NEFD(  1 0) 

CHARACTER*  8  PNAME(7) 

DIMENSION  VAL(7),IPOINT(7),IVARY(7) 

COMMON  /BMOD/XM,V,PC0,SF0,XEFT,NEFD,TBCC,VR,VA,CA 
COMMON  /BMODA/IFLAG,TF 

COMMON  /VALUE/  VAL,IPOINT,IVARY,PNAME,NFIX 
C 

DELTAT=(CTIME-TIME)/ITSTEPS 


c 

C  This  version  of  FDSOLVE  solves  the  compartment 
C  equations  using  backward  Euler  time  differencing 
C  to  avoid  stability  problems.  The  first-order 
C  accuracy  should  be  OK  since  it  matches  the  first-order 
C  spatial  accuracy 

C 

C  Note  a  Gauss-Seidel  solver  is  used  to  solve  the 

C  tridiagonal  matrix  system  that  results  from  the 

C  backward  Euler  application. 

C 

C  This  version  of  fdsolvem  set  up  for: 

C  BATCH  systems 

C  Fitting  multiple  experiments 

C 

C  Set  initial  conditions  from  XIC  input 

C  Also  use  this  for  initial  guess  for  Gauss-Seidel  solver 

C 

RKD=VAL(1) 

DEFF=VAL(2) 

DPATH=VAL(3) 

RLAM=VAL(4) 

XSS0=VAL(5) 

FAST=VAL(6) 

DO  10 1=1  ,NFRAC 
SOLD(I)=XIC(I) 

SNEW(I)=XIC(I) 

10  CONTINUE 

C 

NCOUNT=l 

C 

C  Set  up  constant  used  to  avoid  extra  calculations 

C  see  Dave's  notes 

K(  1  )=2.0*DEFF*DELT  AT/  (DELT  A(  1  )*  *2) 

DO  20  J=2,NFRAC 

K(J)=2.0*DEFF*DELTAT/((DELTA(J- 1  )+DELTA(J))*DELTA(J)) 
20  CONTINUE 

C 

C  Set  up  matrix  constants 

DO  30  J=1,NFRAC-1 
C  Diagonal  term 

A(J)=1 .0+K(J)+FN(J+l)*K(J+l  )/FN(J) 

30  CONTINUE 

C 

A(NFRAC)=1 .0+K(NFRAC) 

DO  40  J=2,NFRAC 
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C  Subdiagonal  term 

B(J)=-K(J)*FN(J)/FN(J- 1 ) 

40  CONTINUE 

C  Note  that  the  super  diagonal  term  is  just 
C  equal  to  -k(j),  j  -  1  ,nfrac- 1 

C 

DO  WHILE  (NCOUNT.LE.ITSTEPS) 

C 

C  Gauss-Seidel  iteration 

IGITER=0 
RESN=1.0 

DO  WHILE((RESN.GT.0.00001).AND.(IGITER.LT.50)) 

C 

SNEW(  1  )=(K(2)*SNEW(2)+SOLD(  1  )+K(  1  )*FN(  1 ) 
*RKD*(1.D0-FAST)*CL)/A(1) 

C 

DO  50  J=2,NFRAC-1 

SNEW(J)=(SOLD(J)+K(J+l)*SNEW(J+l)-B(J)*SNEW(J- 1  ))/A(J) 

C 

50  CONTINUE 

SNEW(NFRAC)=(SOLD(NFRAC)-B(NFRAC)*SNEW(NFRAC- 1  ))/A(NFRAC) 
C 

C  Compute  residual 
C 

RESN=0.0 

RES(  1  )=A(  1  )*SNEW(  1  )-K(2)*SNEW(2)- 

SOLD(  1  )-K(  1  )*FN(  1  )*RKD*CL*(1  .DO-FAST) 

C 

RESN=RESN+(RES(1))**2 

C 

DO  60  J=2,NFRAC-1 
C 

RES(  J)=A(  J)  *  SNEW  ( J)-K(  J+ 1  )*  SNEW  ( J+ 1 )+ 
B(J)*SNEW(J-l)-SOLD(J) 

RESN=RESN+(RES(J)**2) 

60  CONTINUE 
C 

RES(NFRAC)=A(NFRAC)*SNEW(NFRAC)+ 

B(NFR  AC)  *  SNEW  (NFR  AC- 1  )-SOLD(NFRAC) 
RESN=RESN+(RES(NFRAC)**2) 

RESN=SQRT(RESN) 

IGITER=IGITER+ 1 
C 

C  End  Gauss-Seidel  loop 
C 


END  DO 
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c 

C  Add  the  sorbed  concentrations  to  update  cl 
C 

ST=SNEW(1) 

DO  70  J=2,NFRAC 
ST=ST+SNEW(J) 

70  CONTINUE 
C 

C  Update  the  aqueous  concentration  using  the 

C  sorbed  concentration  solutions 

C 

CL=CL+(XM(IEXP)/(VS+RKD*FAST*XM(IEXP)))*(STOLD-ST) 
ANS(NCOUNT)=CL 
ANS2=(ST+FAST*RKD*CL)/CL 
IF  (IFLAG.EQ.O)  ANS(NCOUNT)=ANS2 
C 

DO  80I=1,NFRAC 
SOLD(I)=SNEW  (I) 

XIC(I)=SNEW  (I) 

80  CONTINUE 
C 

STOLD=ST 
NCOUNT=NCOUNT+l 
END  DO 
RETURN 
END 
C 

SUBROUTINE  INTERP(XT,TIME,CTIME,  AN  S  ,ITSTEPS  ,FR) 
DOUBLE  PRECISION  ANS(0:500001),TIME,CTIME,XT,DELTAT 
DOUBLE  PRECISION  TIMEC 
INTEGER  ITSTEPS 
C 

DELTAT=(CTIME-TIME)/ITSTEPS 

NCOUNT=0 

TIMEC=TIME 

DO  WHILE  (XT.GT.TIMEC) 

TIMEC=TIMEC+DELT  AT 
NCOUNT=NCOUNT+l 
END  DO 

C  FR=(ANS(NCOUNT)+ANS(NCOUNT- 1  ))/2.D0 

FR=(ANS(NCOUNT)-ANS(NCOUNT  - 1 ))  *  (XT -TIMEC+DELT  AT) 
./DELTAT+ANS(NCOUNT- 1 ) 

IF  (XT.EQ.TIME)  FR=ANS(0) 

IF  (XT.EQ.CTIME)  FR=ANS(ITSTEPS) 

RETURN 

END 
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c 

SUBROUTINE  INTF(RLAM, DELTA, DMAX,NFRAC,FN) 

DOUBLE  PRECISION  DELTA(501),FN(501),RLAM,DMAX 
INTEGER  NFRAC 
C 

FAC=(RLAM+1  )/(DMAX**(RLAM+l )) 

1=1 

X=0.0 

DO  WHILE  (I.LT.NFRAC) 

C  Apply  Simpson's  rule 

C  FN(I)=DELTA(I)*FAC/6.0*((DMAX-X)**RLAM+ 

C  .4.0*  (DM  AX-(X+DELT  A(I)/2))  *  *RLAM+ 

C  .(DMAX-(X+DELTA(I)))**RLAM) 

FN(I)=(((DMAX-X)**(RLAM+1.D0))-(DMAX-(X+DELTA(I)))** 

.  (RLAM+1  .D0))/(DMAX**(RLAM+1  .DO)) 

X=X+DELTA(I) 

1=1+1 
END  DO 

C  FN(NFRAC)=FAC*DELTA(NFRAC) 

FN(NFRAC)=(DELTA(NFRAC)**(RLAM+1  .DO))/ (DMAX*  *(RLAM+ 1  .DO)) 

RETURN 

END 
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n  n  n  non 


Gamma  Function  Version 


PROGRAM  MODEL 
C 

C  This  program  uses  the  MSS  model  for  sorption  in  a  batch  system 
C  experiencing  liquid  phase  perturbations  and  multiple  experiments. 

C 

C  THIS  VERSION  USES  A  GAMMA  DISTRIBUTION 
C 

C  This  version  calculates  the  sum  of  square  of  error 

C  between  simulated  and  measured  data. 

C 

C  Ed  Heyse/Dave  Coulliette  3  Dec  97 
C  Air  Force  Institute  of  Technology 

C 

REAL  RKD,DEFF,  ALPHA,  BETA,  XSSO,  FAST 

DOUBLE  PRECISION  XM(10),V(10),PC0(10),XEFT(10),SF0(10) 

DOUBLE  PRECISION  TBCC(220),VR(220),VA(220),CA(220) 

DOUBLE  PRECISION  X(250),XALPHA(200),XBETA(200) 
DOUBLE  PRECISION  RFASTV,XMADD,XMREM,XNEWM,TF 
DOUBLE  PRECISION  XIC(501),DELTA(501),  FN(501) 

DOUBLE  PRECISION  COBS,DMAX,COBS2 
DOUBLE  PRECISION  STOLD,VS,CL 

DOUBLE  PRECISION  ANS(0: 10000001), TIME, CTIME,PTIME 
DIMENSION  IX(250),Y  (250),F(250) 

DIMENSION  IFC(220)  ,NEFD(  1 0) 

DIMENSION  V  AL(7),IPOINT(7)  ,1 V  ARY  (7) 

COMMON  /BMOD/XM,V,PC0,SF0,XEFT,NEFD,TBCC,VR,VA,CA 
COMMON  /BMODA/IFLAG,TF 
COMMON  /VALUE/  VAL,IPOINT,IVARY,NFIX 
COMMON  /BDAT/  IX,X,NEXP 
CHARACTER*  80  TIT,LIN1  ,LIN2,LIN3,LIN4 
INTEGER  FDOMAIN 
CHARACTER*  15  INFIL,OUTl 

Initialize  parameters 

COBS=3D0.d0 
NFIX=3D6 

Names  of  input  and  output  files 

WRITE(*,*)  'NAME  OF  INPUT  FILE?' 

READ(5,1 1 1 1)INFIL 
1111  FORMAT(A15) 

WRITE(6,*)'NAME  OF  OUTPUT  FILE?' 


nnn 


READ(5,1 1 1 1)0UT1 

0PEN(UNIT=3  ,FTLE=INFIL,ST  ATUS='OLD') 
OPEN(UNIT=10,FILE=OUT1,STATUS='NEW') 

Read/write  input  data 

READ(3,'(A)')  TIT 
WRITE(  1 0,'(  A)')TIT 
C 

C  Read  parameters  to  be  fitted 
C 

C  RKD=  partition  coefficient,  KP, 

C  DEFF  =  Coefficient  of  Diffusion,  D_eff  cmA2/sec 
C  ALPHA  =  Shape  parameter  for  gamma  distribution 
C  BETA  =  scaling  parameter  for  gamma  distribution 
C  XSSO  =  Initial  Concentration  in  slow  soil  sites,  Sso,  ug/g 

C  FAST  =  Fraction  of  equilibrium  sorption  sites 

C 

READ(3,'(A)')  LIN1 

READ(3,*)  RKD,  DEFF,  ALPHA,  BETA,  XSSO,  FAST 

C  WRITE(10,'(A)')LIN1 

C  WRITE( 10,590)  RKD, DEFF, ALPHA, BETA, XSSO, FAST 

590  FORM  AT(6(  1  X,6E  1 2.5)) 

C 

C  Read  Constants 

C 

C  XM  =  Mass  of  soil  in  reactor,  g 
C  V  =  Volume  of  solvent  in  reactor  at  time  zero,  mL 
C  PCO  =  initial  concentration  in  liquid,  mg/1 

C  SFO  =  experiment  specific  initial  concentration  in  solid 

C 

C  Note:  Intial  cone  in  solid  phase  =  XSS0*SF0(IEXP)  if  initial  sorbed 

C  cone  is  the  same  for  all  experiments,  make  XSSO  =  S(t-O)  and 

C  SFO(IEXP)  =  1 .0  if  initial  sorbed  cone  is  different  for  each  experiment, 

C  make  XSSO  =  1 .0  and  SFO(IEXP)  =  S(t=0,  exp=IEXP) 

C 

C  NEFD  =  Number  of  boundary  condition  changes 
C  IFC  =  experiment  number 

C  TBCC  =  time  for  boundary  condition  to  end  (sec) 

C  VS(J)  =  volume  of  liquid  in  rx  for  BC  period  J,  ml 

C  VR(J)  =  volume  of  liquid  removed  at  end  of  BC  period  J,  ml. 

C  VA(j)  =  volume  of  liquid  added  at  end  of  BC  period  j,  ml. 

C  CA(J)  =  concentration  in  liquid  added  at  end  of  BC  period  J,  mg/1. 

C  XIC  =  sorbed  concentration  in  each  slow  compartment,  |ig/g 
C  CL  =  current  liquid  concentration,  mg/1 
C 
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READ(3,'(A)')  LIN2 
READ(3,*)  (IVARY (I),I=  1  ,NFIX) 

C 

C  READ  IN  HIGH  AND  LOW  values  for  alpha  and  beta 
C 

READ(3,*)NPAIR 
DO  7  IPAIR=1,NPAIR 

READ(3,*)  XALPHA(IPAIR),  XBETA(IPAIR) 

7  CONTINUE 

READ(3,'(A)')  LIN3 

READ(3,*)NEXP 

NABC=0 

C  WRITE(10,'(A)')  LIN2 

C  WRITE(  10,*)  (IVARY(I),I=  1  ,NFIX) 

C  WRITE(10,'(A)')  LIN3 

C  WRITE(10,*)NEXP 

DO  17  I=1,NEXP 

READ(3,*)  XM(I),V(I),PC0(I),SF0(I),NEFD(I) 

C  WRITE(  10,591)  XM(I),  V  (I)  ,PC0(I)  ,SF0(I)  ,NEFD(I) 

NABC=NABC+NEFD(I) 

17  CONTINUE 

591  FORMAT(1X,4E12.5,I10) 

DO  6  I=1,NABC 

READ(3 , *)  IFC(I) ,TBCC(I) , VR(I) , V A(I) ,C A(I) 

C  WRITE(  1 0,592)  IFC(I),TBCC(I),VR(I),VA(I),CA(I) 

6  CONTINUE 

592  FORMAT(  1  X,I5 ,5E1 2.5) 

C 

C  Read  numerical  parameter 

C  IFLAG=0  for  sorbed/liquid  concentration 

C  IFLAG=1  for  liquid  concentration 

C  TF  =  time  factor  =  time  of  time  step/time  for  mass  transfer 
C  Try  TF  =  0.001 

C 

READ(3,*)  IFLAG,TF 
NFRAC1=50 
NFRAC2=100 
C 

C  WRITE(  1 0,*)  IFLAG,TF,NFRAC  1  ,NFRAC2 

CONTINUE 
C 

C  Read  Data  to  be  fitted 

C 

C  X  =  time,  seconds 

C  Y  =  observation  at  time  X,  see  IFLAG. 

C 


no  n  o  o  n  n 


READ(3,'(A)')  LIN4 
C  WRITE(10,'(A)')  LIN4 
NOB=0 
DO  3 1=1,250 

READ(3,*,END=4)  IX(I),X(I),Y(I) 

IF  ((I.GE.2).  AND.(IX(I).NE.IX(I- 1 )))  XEFT(IX(I-1))=X(I-1) 

3  NOB=NOB+l 

4  CONTINUE 
XEFT(IX(NOB))=X(NOB) 

AT  THIS  POINT  WE  CYCLE  THROUGH  VARIOUS  VALUES  OF 

PARAMETERS  TO  CACULATE  SSQ  BETWEEN  OBSERVED 
AND  PREDICTED  RESPONSES 

DO  300  IPAIR2=1  ,NPAIR 
ALPHA=XALPHA(IPAIR2) 

BETA=XBETA(IPAIR2) 

SET  VAL  PARAMETERS  TO  PASS  TO  SUBS 
VAL(1)=RKD 
VAL(3)=ALPHA 
VAL(2)=DEFF 
VAL(4)=BETA 
VAL(5)=XSS0 
VAL(6)=FAST 
C 

CDLC 

C  Begin  new  program  for  batch  solution 

C 

C  NFRAC  is  the  number  of  compartments 
C 

C  NFRAC=(RLAM+2.dO)*  10+1 

C  IF  (NFRAC.LT.50)  NFRAC=50 

C  IF  (NFRAC.GT.500)  WRITE(6,*)  'NFRAC  = NFRAC 

C  IF  (NFRAC.GT.500)  WRITE(10,*)  'NFRAC  = ',  NFRAC 

C 

C  Loop  through  each  experiment  and  set  up  appropriate 

C  initial  conditions 

C 

NSUM=0 

IT=1 

DO  200  IEXP=1,NEXP 
C 

CL=(PC0(IEXP)  *  V(IEXP)+XM(IEXP)  *  XSS0*  SF0(IEXP)*FAST)/ 

.  ( V (IEXP)+RKD*XM(IEXP)  *FAST) 

VS=V(IEXP) 


c 

C  Determine  upper  limit  of  integration  of  gamma  distribution,  DMAX 
C 

SD=-16.661*ALPHA**3.0+38.55*ALPHA**2-31.755*ALPHA+16.066 
IF  (ALPHA.GE.0.59369)SD=- 1 .2844*LOG(ALPHA)+6.6449 
DMAX=BETA*(ALPHA+SD*ALPHA**0.5d0) 

C 

C  Input  the  domain  partitioning.  Delta(I)  is  the  delta  for  each  compartment. 
C 

NFRAC=NFRAC  1 +NFRAC2 

DO  10  I=1,NFRAC1 

DELT  A(I)=ALPH  A*  BET  A/NFR  AC  1 

10  CONTINUE 

DO  11  I=NFRAC1+1,NFRAC 
DELTA(I)=(DMAX-ALPHA*BETA)/NFRAC2 

11  CONTINUE 
C 

C  Subroutine  INTF  generates  the  F(i)'s  by  integrating  f(delta) 

C  using  the  gamma  distribution 
C 

CALL  INTF(ALPHA, DELTA, BETA, NFRAC,FN) 

C 

C  Set  the  initial  condition  within  the  immobile  zone 
C 

C  STOLD  is  the  sum  of  the  S  for  the  diffusion  compartments 

C 

STOLD=0.0 
DO  20  I=1,NFRAC 

XIC(I)=XS  SO*  SFO(IEXP)  *FN(I)*  ( 1  .DO-FAST) 
STOLD=STOLD+XIC(I) 

20  CONTINUE 

C 

C  Loop  through  each  flow  domain 

C 

FDOMAIN=l 

TIME=0.0 

C  IF  (IFLAG.EQ.l)  THEN 

C  WRITE(  1 0,777)IEXP,TIME,PC0(IEXP) 

C  ELSE 

C  WRITE(  1 0,777)IEXP,TIME,XSS0*SF0(IEXP)/PC0(IEXP) 

C  END  IF 

C 

C  IT=1 

C  Note  the  final  time  is  X(NOB) 

DO  400  FDOMAIN=l  ,NEFD(IEXP)+1 


C 
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C  ITSTEPS  is  the  number  of  time  steps  per  flow  domain 

C  CTIME  is  the  cumulative  solution  time 
C 

C  Calculate  ANS  for  beginning  of  boundary  condition 

C 

IF  (IFLAG.NE.O)  ANS(0)=CL 

IF  (IFLAG.EQ.O)  ANS(0)=(STOLD+FAST*RKD*CL)/CL 
C 

CTIME=TBCC(FDOMAIN+NSUM) 

IF  (FDOMAIN.EQ.(NEFD(IEXP)+ 1 ))  CTIME=XEFT(IEXP) 

C 

C  Calculate  retardation  factor  for  slow  sites 
C 

RSLOW=(XM(IEXP)*RKD*(l.dO-FAST))/VS 

C 

C  Calculate  the  number  of  time  steps  for  the  BC 

C 

ITHALF=(CTIME-TIME)*RSLOW*FN(  1 )  *  DEFF/(TF*  DELT  A(  1  )**2.d0)+l 
ITSTEPS=ITHALF*  2 
C 

C  Set  minimum  value  of  ITSTEPS 
C 

IF  ((CTIME-TIME)  .GT.(XEFT (IEXP)/200.d0))  GOTO  301 
GOTO  302 

301  IF  (ITSTEPS.LT.  100)  ITSTEPS=100 

302  CONTINUE 
C 

C  Warn  if  number  of  time  steps  exceeds  array  size 
C 

IF  (ITSTEPS .GT.  10000000)  WRITE(6,*)  'ITSTEPS  = ITSTEPS 
IF  (ITSTEPS.GT.  10000000)  WRITE(10,*)  'ITSTEPS  = ',  ITSTEPS 
IF  (ITSTEPS.GT.  10000000)  GOTO  2222 
C 

C  FDSOLVE  computes  the  soln  from  time  to  ctime  for  itsteps  and 
C  passes  the  soln  back  as  array  ANS(ITSTEPS) 

C 

CALL  FDSOLVE(TIME, CTIME, ITSTEPS, ANS, XIC, 

.  FN,DELTA,CL,STOLD,VS,IEXP,NFRAC) 

C 

C  Write  prediction  data  to  output  file 
C 

PTIME=TIME 

DELTAT=(CTIME-TIME)/ITSTEPS 
C  DO  7711=0,10 

C  CALCT=PTIME 

C  IF  (II.EQ.0)  CALCT =PTIME+ 1  .d0 
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C  IF  (II.EQ.  10)  CALCT=PTIME-l.dO 

C  CALL  INTERP(CALCT, TIME, CTIME,ANS,ITSTEPS, COBS) 

C  WRITE( 1 0,777)IEXP,PTIME,COBS,ITSTEPS 

C  WRITE(  1 0,777)IEXP,C  ALCT,COB  S 

C  PTIME=PTIME+0. 1  dO*(CTIME-TIME) 

777  FORMAT  ( 1  X,I6,2(E20.9)) 

C  77  CONTINUE 

C 

C  Calculate  predicted  values  corresponding  to  observed  times 

C 

DO  WHILE  ((X(IT).LE.CTIME).AND.(X(IT).GE.TIME) 
.AND.(IX(IT).EQ.IEXP)) 

CALL  INTERP(X(IT), TIME, CTIME,ANS,ITSTEPS,COBS2) 
F(IT)=COBS2 
IT=IT+1 
END  DO 
C 

TIME=CTIME 

C 

C  Adjust  aqueous  volume  and  new  aqueous  cone. 

C 

RFASTV=XM(IEXP)*RKD*FAST+VS 

XMADD=VA(FDOMAIN+NSUM)*CA(FDOMAIN+NSUM) 

XMREM=VR(FDOMAIN+NSUM)*CL 

XNEWM=RFASTV*CL+XMADD-XMREM 

CL=XNEWM/(RFASTV- 

VR(FDOMAIN+NSUM)+VA(FDOMAIN+NSUM)) 

VS=VS-VR(FDOMAIN+NSUM)+VA(FDOMAIN+NSUM) 

C 

C  Statement  400  =  end  of  do  loop  for  BC  change 
C 

400  CONTINUE 

C 

C  Statement  200  =  end  of  do  loop  for  experiment 
C 

NSUM=NSUM+NEFD(IEXP) 

200  CONTINUE 
C 

SSQ=0. 

DO  253  IS=l,NOB 

C  WRITE(  10,777)  IX(IS),X(IS),Y(IS) 

SSQ=SSQ+(Y(IS)-F(IS))**2. 

253  CONTINUE 

WRITE(  10,254)  ALPHA, BETA, SSQ 
WRITE(6,254)  ALPHA, BETA, SSQ 

254  FORMAT( IX, 'Alpha  =  ',E16.5,'  Beta  =  ',E16.5,'  SSQ  =  ',E16.5) 


300  CONTINUE 
C 

2222  CONTINUE 
STOP 
END 
C 

SUBROUTINE  INTF(ALPHA,  DELTA,  BETA,  NFRAC,FN) 
DOUBLE  PRECISION  DELTA(501),FN(501),SUM 
REAL  ALPHA, BETA 
INTEGER  NFRAC 
EXTERNAL  FUNC 
C 

C  SUBROUTINE  INTF  INTEGRATES  DISTRIBUTION  FUNCTION 

Q  *************Q^MMA************** 

C  FOR  EACH  COMPARTMENT  I 

C 

SUM=0.d0 

1=1 

C 

C  Lower  Limit  of  integration: 

C  alpha  >=  1 :  LL  =  0 

C  alpha  <  1 :  Start  from  compartment  2,  est 

C  compartment  1  from  0.99955  -  sum(F2..FN) 

C 

IF  (ALPHA.LT.l  .0)X=DELTA(1) 

IF  (ALPHA.LT.l. 0)1=2 
IF  (ALPHA.LT.  1.0)GOTO  10 
IF  (ALPHA.GE.  1 .0)X=0.d0 

C  IF  (ALPHA.LT.  1 .0)X=ALPHA*BETA*  1 0.**(-3.8/ ALPHA-0.2) 

C  X=(l  .OE-08)*DELTA(1) 

10  CONTINUE 

DO  WHILE  (I.LE.NFRAC) 

C 

X 1  =X+DELT  A(I) 

C  INTEGRATE  FROM  X  TO  X 1  USING  SUBROUTINE  QROMB 

C 

CALL  QROMB(FUNC, X, XI , ALPHA, BETA, Tl) 

FN(I)=T1 

C 

SUM=SUM+FN(I) 

1=1+1 
X=X1 
END  DO 

IF  (ALPHA.LT.  1.0)  FN(1)=0.99955-SUM 

RETURN 

END 


c 

FUNCTION  FUNC(X,ALPHA,BETA) 

REAL  ALPHA, BETA, X 
C 

T1  =EXP(GAMMLN(ALPHA)) 

FUNC=(BETA**(-ALPHA))*(X**(ALPHA-1))*(EXP(-X/BETA))/T1 

RETURN 

END 

C 

SUBROUTINE  qromb(func,  a,  b,  alpha,  beta,  ss) 

INTEGER  JMAX,JMAXP,K,KM 
REAL  a,  b,  func,  ss,  EPS 
EXTERNAL  func 

PARAMETER  (EPS=l.e-7,  JMAX=30,  JMAXP=JMAX+1,  K=8,  KM=K-1) 
CU  USES  polint,trapzd 
INTEGER  j 

REAL  dss,  h(JMAXP),  s(JMAXP) 
h(l)=L 

do  11  j=l,JMAX 

call  trapzd(func,  a,  b,  alpha,  beta,  s(j),  j) 
if  (j.ge.K)  then 

call  polint(h(j-KM),s(j-KM),K,0.,ss,dss) 
if  (abs(dss).le.EPS*abs(ss))  return 
endif 

s(j+l)=s(j) 

h(j+l)=0-25*h(j) 

1 1  continue 

pause  'too  many  steps  in  qromb' 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  $!6)$#15[)Nu. 

C 

SUBROUTINE  polint(xa,ya,n,x,y,dy) 

INTEGER  n,NMAX 
REAL  dy,x,y,xa(n),ya(n) 

PARAMETER  (NMAX=10) 

INTEGER  i,m,ns 

REALden,dif,dift,ho,hp,w,c(NMAX),d(NMAX) 

ns=l 

dif=abs(x-xa(l)) 

do  1 1  i=l,n 

dift=abs(x-xa(i)) 

if  (dift.lt.dif)  then 

ns=i 

dif=dift 

endif 

c(i)=ya(i) 
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u  u 


d(i)=ya(i) 

1 1  continue 
y=ya(ns) 
ns=ns-l 

do  13  m-l,n-l 

do  12  i=l,n-m 

ho=xa(i)-x 

hp=xa(i+m)-x 

w=c(i+l)-d(i) 

den=ho-hp 

if(den.eq.O.)pause  'failure  in  polint' 

den=w/den 

d(i)=hp*den 

c(i)=ho*den 

12  continue 

if  (2*ns.lt.n-m)then 

dy=c(ns+l) 

else 

dy=d(ns) 

ns=ns-l 

endif 

y=y+dy 

13  continue 
return 
END 

(C)  Copr.  1986-92  Numerical  Recipes  Software  $!6)$#15[)Nu. 

SUBROUTINE  trapzd(func, a, b, alpha, beta, s,n) 

INTEGER  n 

REAL  a,b,s,func,alpha,beta 
EXTERNAL  func 
INTEGER  it,j 
REAL  del,sum,tnm,x 
if  (n.eq.l)  then 

s=0.5*(b-a)*(func(a, alpha, beta)+func(b, alpha, beta)) 
else 

it=2**(n-2) 
tnm=it 

del=(b-a)/tnm 
x=a+0.5*del 
sum=0. 
do  11  j=l,it 

sum=sum+func(x, alpha, beta) 
x=x+del 
1 1  continue 

s=0.5*(s+(b-a)*sum/tnm) 
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no  no 


endif 

return 

END 

C  (C)  Copr.  1986-92  Numerical  Recipes  Software  $!6)$#15[)Nu. 

FUNCTION  gammln(xx) 

REAL  gammln,xx 
INTEGER j 

DOUBLE  PRECISION  ser,stp,tmp,x,y,cof(6) 

SAVE  cof,stp 

DATA  cof  ,stp/76. 1 8009 1 72947 1 46d0, -86.50532032941 677d0, 

*24.0 1 40982408309 1  d0,- 1 .23 1 739572450 1 55d0,.  1 208650973866 1 79d-2, 

*-.5395239384953d-5,2.5066282746310005dO/ 

x=xx 

y=x 

tmp=x+5.5dO 

tmp=(x+0.5d0)*log(tmp)-tmp 
ser=l  .0000000001 9001 5d0 
do  11  j— 1,6 
y=y+l  .d0 
ser=ser+cof(j)/y 
1 1  continue 

gammln=tmp+log(stp*ser/x) 

return 

END 

(C)  Copr.  1986-92  Numerical  Recipes  Software  $!6)$#15[)Nu. 

SUBROUTINE  FDSOLVE(TIME,CTIME,ITSTEPS,ANS,XIC, 

.FN, DELTA, CL, STOLD,VS,IEXP,NFRAC) 

IMPLICIT  NONE 

DOUBLE  PRECISION  RKD,DEFF, ALPHA, BETA, XSS0, FAST, TF 
DOUBLE  PRECISION  XM(10),V(10),PC0(10),SF0(10),XEFT(10) 
DOUBLE  PRECISION  TBCC(220),VR(220),VA(220),CA(220) 
DOUBLE  PRECISION  ANS(0:10000001),TIME,CTIME,STOLD 
DOUBLE  PRECISION  XIC(501),DELTA(501),FN(501),CL,VS 
DOUBLE  PRECISION  K(501),RES(501),A(501),B(501),RESN,ANS2 
INTEGER  ITSTEPS,NFRAC 

DOUBLE  PRECISION  SNEW(501),SOLD(501),ST,DELTAT 
DIMENSION  NEFD(10) 

DIMENSION  V  AL(7)  ,IPOINT  (7)  ,1 V  ARY  (7) 

COMMON  /BMOD/XM,V,PC0,SF0,XEFT,NEFD,TBCC,VR,VA,CA 
COMMON  /BMODA/IFLAG,TF 
COMMON  /VALUE/  VAL,IPOINT,IVARY,NFIX 
C 

DELTAT=(CTIME-TIME)/ITSTEPS 
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ng  o  o  g  n  on 


c 

C  This  version  of  FDSOLVE  solves  the  compartment  equations  using  backward 
C  Euler  time  differencing  to  avoid  stability  problems.  The  first-order  accuracy 
C  should  be  OK  since  it  matches  the  first-order  spatial  accuracy. 

C 

C  Note  a  Gauss-Seidel  solver  is  used  to  solve  the  tridiagonal  matrix  system 
C  that  results  from  the  backward  Euler  application 
C 

C  This  version  of  fdsolvem  set  up  for: 

C  BATCH  systems  fitting  multiple  experiments 

C 

C  Set  initial  conditions  from  XIC  input 

C  Also  use  this  for  initial  guess  for  Gauss-Seidel  solver 

C 

RKD=VAL(1) 

DEFF=VAL(2) 

ALPHA=VAL(3) 

BETA=VAL(4) 

XSS0=VAL(5) 

FAST=VAL(6) 

DO  10I=1,NFRAC 
SOLD(I)=XIC(I) 

SNEW(I)=XIC(I) 

10  CONTINUE 

C 

NCOUNT=l 

C 

C  Set  up  constant  used  to  avoid  extra  calculations 

C  see  Dave's  notes 

K(  1  )=2.0*  DEFF*  DELT  AT/(DELTA(  1 )  *  *  2) 

DO  20  J=2,NFRAC 

K(J)=2.0*DEFF*DELTAT/((DELTA(J- 1  )+DELT  A( J))  *  DELT  A(  J)) 

20  CONTINUE 

Set  up  matrix  constants 
DO  30  J=1,NFRAC-1 
Diagonal  term 

A(J)=1 .0+K(J)+FN(J+l  )*K(J+1  )/FN(J) 

CONTINUE 

A(NFRAC)=1 .0+K(NFRAC) 

DO  40  J=2,NFRAC 
Subdiagonal  term 
B(J)=-K(J)*FN(  J)/FN(J- 1 ) 

CONTINUE 

Note  that  the  super  diagonal  term  is  just  equal  to  -k(j),  j=l  ,nffac-l 


68 


c 

DO  WHILE  (NCOUNT.LE.ITSTEPS) 

C 

C  Gauss-Seidel  iteration 
IGITER=0 
RESN=1.0 

DO  WHILE((RESN.GT.0.00001).AND.(IGITER.LT.50)) 

C 

SNEW(l)=(K(2)*SNEW(2)+SOLD(l)+K(l)*FN(l) 

*RKD*(  1  .D0-FAST)*CL)/  A(  1 ) 

C 

DO  50  J=2, NFRAC- 1 

SNEW(J)=(SOLD(J)+K(J+l)*SNEW(J+l)-B(J)*SNEW(J-l))/A(J) 

C 

50  CONTINUE 

SNEW(NFRAC)=(SOLD(NFRAC)-B(NFRAC)*SNEW(NFRAC- 1  ))/A(NFRAC) 
C 

C  Compute  residual 
C 

RESN=0.0 

RES(  1  )=A(  1  )*SNEW  ( 1  )-K(2)*SNEW  (2)- 

SOLD(  1  )-K(  1 )  *FN(  1 )  *RKD*  CL*  ( 1  .DO-FAST) 

C 

RESN=RESN+(RES(1))**2 

C 

DO  60  J=2,NFRAC-1 
C 

RES(J)=A(J)*SNEW(J)-K(J+1)*SNEW(J+1)+ 

B(J)*SNEW(J-l)-SOLD(J) 

RESN=RESN+(RES(J)*  *2) 

60  CONTINUE 
C 

RES(NFRAC)=A(NFRAC)*SNEW(NFRAC)+ 

B  (NFR  AC)  *  SNEW  (NFRAC- 1  )-SOLD(NFRAC) 
RESN=RESN+(RES(NFRAC)**2) 

RESN=SQRT(RESN) 

IGITER=IGITER+ 1 
C 

C  End  Gauss-Seidel  loop 
C 

END  DO 
C 

C  Add  the  sorbed  concentrations  to  update  cl 
C 

ST=SNEW(1) 

DO  70  J=2, NFRAC 
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ST=ST+SNEW(J) 

70  CONTINUE 
C 

C  Update  the  aqueous  concentration  using  the 

C  sorbed  concentration  solutions 

C 

CL=CL+(XM(IEXP)/(VS+RKD*FAST*XM(IEXP)))*(STOLD-ST) 
ANS(NCOUNT)=CL 
ANS2=(ST+FAST*RKD*CL)/CL 
IF  (IFLAG.EQ.O)  ANS(NCOUNT)=ANS2 
C 

DO  80  I=1,NFRAC 
SOLD(I)=SNEW(I) 

XIC(I)=SNEW  (I) 

80  CONTINUE 
C 

STOLD=ST 

NCOUNT=NCOUNT+l 

END  DO 

RETURN 

END 

SUBROUTINE  INTERP(XT, TIME, CTIME,ANS,ITSTEPS,FR) 

DOUBLE  PRECISION  ANS(0: 10000001), TIME, CTIME,XT,DELTAT,FR 
DOUBLE  PRECISION  TIMEC 
INTEGER  ITSTEPS 
C 

DELTAT=(CTIME-TIME)/ITSTEPS 

NCOUNT=0 

TIMEC=TIME 

DO  WHILE  (XT.GT.TIMEC) 

TIMEC=TIMEC+DELT  AT 
NCOUNT=NCOUNT+l 
END  DO 

C  FR=(ANS(NCOUNT)+ANS(NCOUNT- 1  ))/2.D0 

FR=(ANS(NCOUNT)-ANS(NCOUNT-l))*(XT-TIMEC+DELTAT) 
./DELTAT+ANS(NCOUNT- 1 ) 

IF  (XT.EQ.TIME)  FR=ANS(0) 

IF  (XT.EQ.CTIME)  FR=ANS(ITSTEPS) 

RETURN 

END 
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APPENDIX  D:  Details  of  RSM  Procedure  for  Radial  Geometry  Function 

This  appendix  details  the  RSM  procedure  followed  for  the  radial  geometry  function 
version  of  the  MSS  model.  Mika’s  (1997)  even  data  set,  condition  5  (Appendix  B)  is  used  in 
the  MSS  model  runs  to  obtain  SSE  values  that  form  the  response  surface. 

The  first  minimum  bias  design  of  experiments  for  the  radial  geometry  function  is 
given  in  Table  D-l.  These  variables  are  used  in  the  statistical  package  SASJMP ™  to  obtain 
the  regression  function,  y  =  0.06247  -  0.007094  -  0.0461 5 -R,  which  is  the  first-order 

approximation  of  the  response,  SSE. 

Table  D-l,  First  22  Factorial  Design  for  Radial  Geometry  Function 


Run  # 

Coded  Variables 

Uncoded  Variables 

Response 

SSE 

X, 

x2 

X 

R 

1 

-VV3 

-VP 

1.2887 

0.1683 

0.01515 

2 

-VP 

VP 

1.2887 

0.2683 

0.07136 

3 

VP 

-VP 

2.2887 

0.1683 

0.06888 

4 

VP 

VP 

2.2887 

0.2683 

0.00344 

The  next  step  is  to  conduct  the  gradient  search  to  determine  possible  directions  of 
improvement.  Since  an  improvement  in  this  case  is  defined  to  be  a  lower  SSE,  y  is 
multiplied  by  (-1)  before  taking  the  partial  derivative  to  obtain  the  gradient  vector.  After  the 
gradient  vector  is  divided  by  its  magnitude,  the  unit  gradient  vector  is  obtained,  which  is  also 
the  step  size  to  take  for  new  experiments  (in  coded  units).  By  use  of  equation  21 ,  these 
coded  units  are  transformed  to  the  original  variables  which  are  run  in  the  MSS  model.  The 
results  are  shown  in  Table  D-2.  Notice  that  the  SSE  increases  before  and  after  point  4, 
indicating  that  this  may  be  the  minimum  response. 
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To  obtain  a  better  first-order  fit,  point  4  is  used  as  the  center  point  for  another  22 
factorial  design,  followed  by  a  gradient  search.  Table  D-3  contains  the  22  factorial  data  as 
well  as  the  gradient  search  results.  Again,  notice  that  point  4  has  the  minimum  response. 

Using  the  data  in  Table  D-3,  an  augmented  first-order  design  is  completed  to 
determine  if  the  fit  will  be  satisfactory.  The  regression  equation  is 

y  =  -0.4195  +0.2585-X.  +  2.1298R-  1.2165A-R 

with  a  correlation  coefficient  (R2)  of  0.7945.  This  coefficient  is  probably  low  due  to 
curvature  of  the  surface.  Since  a  higher  R2  is  desired,  a  second-order  response  will  be  fit. 


Table  D-2,  Gradient  Search  Results 

Gradient  Vector,  b  -0  022898  0.751635 

Magnitude  of  6  0.7519837 

Unit  Gradient  Vector  -0.03045026  0.99953628 


Coded  Variables  | 

Uncoded  Variables  ] 

Response 

Run# 

*1 

*2 

X 

R 

SSE 

Center  Point 

0 

0 

1.5 

0.175 

0.01930020981 

Step  Size 
Point  1 

-0.03045026 

-0.57735027 

0.99953628 

-0.57735027 

-0.01522513 

1.21132487 

0.07496522 

0.13169873 

0.05666088122 

Point  2 

-0.57735027 

0.57735027 

1.21132487 

0.21830127 

0.01561637710 

Point  3 

0.57735027 

-0.57735027 

1.78867513 

0.13169873 

0.09393009159 

Point  4 

0.57735027 

0.57735027 

1.78867513 

0.21830127 

0.00478763914 

Point  5 

-0.03045026 

0.99953628 

1.48477487 

0.24996522 

0.02344801782 

Point  6 

-0.06090052 

1.99907257 

1.46954974 

0.32493044 

0.13963631887 

Point  7 

-0.09135078 

2.99860885 

1.45432461 

0.39989566 

0.30743289973 

Point  8 

-0.12180104 

3.99814513 

1.43909948 

0.47486088 

0.49175314901 

A  central  composite  design  (CCD)  is  constructed  using  point  4  from  Table  D-3  as  the 
center  point  with  one  center  run.  Table  D-4  contains  the  data  as  input  to  the  SAS  JMP ™ 
statistical  package.  With  this  data,  a  response  surface  model  is  constructed,  which  results  in 
a  second-degree  polynomial  approximation  of  the  response  surface.  The  summary  of  fit 
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information  for  this  model  is  presented  in  Table  D-5.  The  correlation  coefficient  indicates  a 
good  fit  to  the  data.  This  fit  is  for  a  small  range  of  X  and  R,  as  can  be  seen  in  Table  D-4. 


Table  D-3,  Second  22  Factorial  Design  and  Gradient  Search 


Gradient  Vector,  b  0.007094  0.046149 

Magnitude  of  6  0.04669106 
Unit  Gradient  Vector  0.15193487  0.98839051 


Coded  Variables 


Uncoded  Variables 

X  _ R 


Run# 

Center  Point 
Step  Size 
Point  1 
Point  2 
Point  3 
Point  4 
Point  5 
Point  6 
Point  7 
Point  8 


*i 

0 

0.15193487 

-0.57735027 

-0.57735027 

0.57735027 

0.57735027 

0.15193487 

0.30386974 

0.45580461 

0.60773948 


X2 

0 

0.98839051 

-0.57735027 

0.57735027 

-0.57735027 

0.57735027 

0.98839051 

1.97678102 

2.96517152 

3.95356203 


1.78867513 

0.13157945 

1.28867513 

1.28867513 

2.28867513 
2.28867513 
1.92025459 
2.05183404 
2.18341349 
2.31499294 


0.21830127 

0.08559713 

0.16830127 

0.26830127 

0.16830127 

0.26830127 

0.3038984 

0.38949553 

0.47509266 

0.56068978 


Response 

SSE 

0.0047876391 

0.0151496775 

0.0713609895 

0.0688815446 

0.0034405227 

0.0388629926 

0.1343904003 

0.2486271589 

0.3638264837 


Table  D-4,  CCD  for  X,  R  Response  Surface 


Run  # 

Coded  Variables 

Uncoded  Variables 

Response, 

SSE 

X, 

x2 

X 

R 

1 

-VP 

-VP 

1.6763 

0.2071 

0.006034 

2 

-VP 

VP 

1.6763 

0.3295 

0.105875 

3 

VP 

-VP 

2.9010 

0.2071 

0.053957 

4 

VP 

VP 

2.9010 

0.3295 

0.006333 

5 

-1 

0 

1.4226 

0.2683 

0.051296 

6 

1 

0 

3.1547 

0.2683 

0.016846 

7 

0 

-1 

2.2887 

0.1817 

0.052650 

8 

0 

1 

2.2887 

0.3549 

0.056158 

9 

0 

0 

2.2887 

0.2683 

0.003441 

Table  D-5,  Summary  Fit  data  of  CCD  for  X,  R  Response  Surface 

R-Square  0.96936 

Adjusted  R-Square  0.91 829 

Root  Mean  Square  Error  0.00971 
Mean  of  Response  0.039 1 7 

Observations  9 
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The  fitted  model  is  now 


y=  0.10856  +  0.059974  -  1.2371-R  +  0.040042  -  0.9831  A-R  +  6.71595R2 

in  the  original  uncoded  units  of  X  and  R.  These  estimates,  the  standard  errors,  and  p-values 
are  presented  in  Table  D-6.  The  p-values  indicate  that  the  only  significant  terms  in  the 
model  are  X2,  XR,  and  R2.  The  effect  of  removing  the  other  terms  is  to  cause  the  correlation 
coefficient  (R2)  to  decrease.  Since  these  terms  are  only  used  to  define  the  response  surface, 


they  will  be  left  in  the  regression  model. 


Table  D-6,  X ,  R  Model  Summary  Statistics 


Term 

Estimate 

Std  Error 

t  Ratio 

p-value 

Intercept 

0.108556 

0.181727 

0.60 

0.5924 

X 

0.059967 

0.077913 

0.77 

0.4975 

R 

-1.2371 

0.868915 

-1.42 

0.2497 

X2 

0.004005 

0.015186 

2.64 

0.0778 

XR 

-0.983094 

0.129505 

-7.59 

0.0047 

R2 

6.715951 

1.518579 

4.42 

0.0215 

A  surface  plot  of  this  second-order  equation  can  be  found  in  Figure  1 1 .  Note  how 
similar  the  response  function  plot  is  to  a  plot  of  actual  model  data  over  the  same  region, 
presented  in  Figure  12. 

The  RSM  procedure  as  outlined  above  is  performed  with  the  gamma  MSS  model  for 
both  the  even  mix  data  set  and  the  Borden  soil  data  set. 
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